Point cloud based tool path generation for corrective machining in ultra-precision diamond turning

The increasing demand for machining non-rotational optical surfaces requires capable and flexible cutting tool path generation methods for ultra-precision diamond turning. Furthermore, the recent interest in on-machine metrology and corrective machining requires efficient as well as accurate algorithms capable to handle point cloud based surface data. In the present work, a new computation method for the tool path generation is proposed that focuses on three-axes corrective machining. It is based on the principle of defining the surface to be machined by a point cloud of given density, since surface measurement data is usually available as point cloud. Numeric approximation techniques are used to compute the surface normal vectors and calculate the resulting positions of the cutting tool path preserving a uniform radial axis motion for face turning. Investigations are performed in order to quantify the error between the calculated tool path and the exact analytical solution. The error dependencies are analyzed regarding the local surface slope and numerical parameters. Error values below 1 nm are achieved. In addition, form deviation results prove the method’s capability for corrective diamond turn machining.


Introduction
Diamond turning and ultra-precision machining technologies today in general offer the capabilities to machine surfaces with form deviations from the nominal surface in the range of 100 nm and roughness values below 5 nm. Besides the use of ultra-precision machine tools, whose structure and components offer high level movements in terms of both precision and accuracy, cutting tools made of single crystal diamonds are employed to achieve the mentioned surface quality. Hence, ultra-precision machining is often referred to as diamond machining [1].
Recently, the popularity of three-axes diamond turning increased in machining of non-rotational symmetric geometries including complex aspheres, lens arrays and freeforms in general. Freeform optics enable new possibilities for designers of optical systems for example in terms of a reduced number of optical elements and an increased optical performance [2]. This leads to a growing demand for precision freeform surfaces and thus forces diamond machining technologies to comply with this request. Corrective machining in combination with on-machine surface metrology is a promising approach for increasing diamond turning capabilities [3]. In view of these developments effective tool path generation is of high importance, especially for corrective machining.
Within the tool path generation the tool geometry, consisting of the nose radius and the rake angle, as well as the machining parameters like feed rates and NC point resolution need to be taken into consideration. Figure 1 shows the generalized steps of machining an optical surface and the accuracy losses throughout the process. Regarding the tool path generation, the aim is to minimize its related accuracy loss in order to decrease the overall accuracy loss along all processing steps. For different surface definition methods or surface data available as measured point cloud in case of corrective machining additional format transformations, e.g., fitting operations, may be necessary.
For the common use case of machining surfaces for optical applications by face turning, the tool path topography mainly follows a spiral. Regarding the tool nose radius compensation, Gong et al. [4] distinguish and discuss two different ways of projecting the spiral tool path onto the surface to be machined. The first is called oscillating-X method, since it can cause oscillating movements in X-direction having the potential to negatively affect the machining process due to high dynamic movements. The second leads to a constant feed in the radial direction and is named as steady-X method. Since constant or slowly changing feed rates prevent the process from any high dynamic movements, the steady-X method is generally preferred, but also more complex in terms of the tool-surface contact points calculation based on a predefined spiral tool path. Several publications discuss the issues in this context and present diverse approaches for the tool path generation for freeform machining. Gong et al. [4] give a comprehensive review of tool path generation techniques in ultra-precision diamond turning. Besides, Brinksmeier et al. [5] focus on tool path generation strategies for freeform surfaces.
A widely used method to describe freeforms is using nonuniform rational B-splines (NURBS), as described in [6,7]. NURBS offer a framework to represent freeforms in a standardized way and to get an analytical description between given sampling points. However, Gong et al. [4] point out that there is no obvious correspondence between the parameter of the NURBS and the coordinate parameters and that analytic parameter conversion is often impossible. Moreover, also stated in [8], the NURBS approximation may not reach the desired accuracy. To overcome the problem it is possible to sample the surface and use a function fitting to obtain an analytical formula, see Scheiding et al. [9]. This process may suffer from accuracy and time consumption, especially in case of sampling noise or uneven point distribution, see [10] and references therein. In [4] it is noted that a similar algorithm is used in NANOCAM 2.0, but which may lead to accuracy losses.
To overcome the drawback of accuracy losses when approximating an analytic surface description by NURBS, Gong et al. [8] present a tool path generation method using symbolic computation for the calculation of the surface normals. However, the finding of the tool-surface contact point still requires numerical computation, since this usually cannot be done analytically for a given spiral path [8].
Moreover, there are machining applications which are forced to handle the surface to be machined defined by(nonuniform or uniform) point clouds. An example is machining related to reverse engineering, where surface information typically consists of coordinate points gained by 3D scanning (as stated in [11,12]). Masood et al. [11] present a direct machining concept for three-axes milling using ballend tools, which avoids the traditional surface fitting process in reverse engineering applications. The point cloud is divided into fields which correspond to one of the parallel oriented milling paths. Each path is then obtained by fitting a curve to the tool center points of the respective field, which are obtained using the normals of the triangulated surface data. The accuracy of the employed 3D scanning equipment is stated to be ±, 0.018 mm. However, no analysis is presented regarding the tool path accuracy compared to the point cloud accuracy and to the method using surface fitting processes.
Chui et al. [12] also present a direct machining approach tailored to milling using ball-end tools. It focuses on the applicability to both regularly and irregularly distributed point clouds. After subdividing the points into fields, whose width depends on the step-over of two adjacent milling paths, the tool radius center points are obtained by geometric calculations determining the contact points between a surface point and the tool. A discussion of the resulting tool path accuracy is not part of their work.
Corrective machining of surfaces is an application, where the actual measured surface form is mainly available as point cloud. In this context, Xiao et al. [13] present an algorithm providing direct tool path generation for the repair machining of aero-engine blades based on a measured point cloud. The approach avoids surface fitting operations as well as the usage of CAM software and is designed to feature five-axes milling. Hence, it is concluded to be more efficient than traditional methods. The accuracy of the employed scanning measurement is able to reach a value of 30 . Besides, algorithms are presented for the region separation of the blades and recognition of damaged areas. The tool path accuracy is not explicitly analyzed, but is stated to be mainly influenced by the point cloud's accuracy and density. As mentioned before, corrective machining is also becoming more popular in ultra-precision diamond turning and is strongly encouraged by the recent interest in on-machine surface metrology [3]. Although the use of a point cloud for the tool path generation is not further analyzed in [5], it is listed in the given overview as a possible strategy for ultraprecision machining. The named fact that numeric computations cannot be avoided entirely for a given spiral path (see [8]) and point cloud form measurement data motivate a point cloud based surface representation throughout the tool path generation process, which, in contrast to a hybrid process according to [8] using both analytic and numeric computations, for example also allows to handle freeform point clouds exported from a CAD design software. These aspects together with the characteristics and benefits of the above-described applications of direct machining presented in [11][12][13] lead to the demand for a tool path generation method combining efficiency in processing point clouds and enhanced accuracy to the needs of diamond turning. Therefore, this work presents a novel simplified method that uses a point cloud and numerical computations to calculate the tool path for 3D diamond turning. It follows the explained and preferred steady-X method to maintain a uniform feed rate in radial direction, is capable to handle various tool rake angles and offers the option to compensate the point cloud by measured form errors in case of corrective machining. Since this strategy implies a certain accuracy loss due to no gapless analytical surface description and the employment of numeric computation methods, a comprehensive accuracy loss analysis is presented as well. The analysis and the method in general follow the aim of calculating tool path points of highest accuracy. This is required for any tool path point in diamond machining, regardless of the point density with which a tool path is programmed.
As part of the analysis every calculation step of the presented approach accompanied by a numerical approximation is investigated regarding its accuracy. Thus, the present work supplements the important analysis, which compares the point cloud approach with the exact solution and hence mitigates the shortcomings in the analysis of [11][12][13]. In contrast to these, the selection of the computational methods and its parameters is deduced and justified by analyzing various conditions regarding the accuracy demands in diamond turning.
Beyond, exemplary form error results prove the usability for freeform corrective machining in combination with the use of on-machine metrology.

Tool path generation
The configuration of the ultra-precision diamond turning machine used for this work is depicted in top view in Fig. 2.
In contrast to the more usual configuration, the spindle is mounted on the Z-axis slide. The two horizontal linear axes are structurally separated and carry either the diamond tool or the workpiece spindle (dependent on the individual machine tool type), which operated as a controlled C-axis enables three-axes machining. Shown is also the stationary Cartesian machine tool coordinate system, in which the X-axis corresponds to the radial linear axis of the machine tool. It serves as the coordinate system used throughout the present work.
As a speciality of diamond turning machines, a movement in Z-direction can be performed by two different modes, namely either by the Z-axis itself (referred to as slow slide servo machining), or by a fast tool servo device (referred to as fast tool servo machining). In this context it is noted that the tool path generation generally is independent of the Z-axis mode used for machining (slow slide or fast tool servo), because both modes move according to the Z-coordinates calculated with respect to the machining coordinate system. The selection of the Z-axis mode takes place when programming the tool path to the machine tool control using the control-specific syntax, but which is not part of the present work.
The machining kinematics in face turning and the typical cutting tool geometry of a circular arc with certain radius cause the contact point between the workpiece surface and the tool to move along the tool nose arc when machining curved geometries. This behavior requires to compensate the nominal spiral tool path in function of the tool nose radius and the surface normal vector at each tool path point. For that, the ability to calculate the surface normal vector at Fig. 2 Axes configuration of the used three-axes diamond turning machine consisting of two horizontal linear axes (X-and Z-axis) and one horizontal rotational axis (C-axis) with stationary right-handed Cartesian machine tool coordinate system. The configuration corresponds to a typical one used in diamond turning each point of the surface to be machined is crucial. Figure 3 illustrates the relation between a surface profile f and the associated tool path.
When calculating the surface normals by analytical methods, additional precomputation for each specific surface equation is necessary. This is avoided by representing the nominal surface by a point cloud as described in Sect. 1. The point cloud represents a flexible basis for the overall tool path generation and can be generated out of any nominal surface description by an initial preprocessing step. Two further steps complete the path generation, namely the calculation of the tool center surface and of the tool path points.

Calculation of surface normals
By machining a general freeform surface using a circular radius cutting tool the tool's radius center point has to move on a surface that has a constant distance to the surface to be machined. Hence, the first step of the tool path generation is the calculation of the normal vectors of the surface to be machined. The surface is considered as a point cloud P consisting of k points w w w i with the coordinates x i , y i and z i according to the definition given by Eq. (1).
As already mentioned, the point cloud can be derived by an initial preprocessing step out of any surface description method. Its properties and the associated effects are discussed in Sect. 3 subsequently.
For the common case of an analytical surface description z = f (x, y) , the surface normal vectors n n n(x, y) can be derived using the known analytical formulas. Although this description method is common for example for aspherical surfaces, it cannot be used for arbitrary freeform surfaces or for the correction of surfaces with measured points, which are usually available as point cloud. Using a point cloud as surface description method overcomes this drawback and represents the basis for a generalized surface normal calculation. Therefore, in a first step, the method described by Hoppe et al. [14] is utilized, where the normal vector to a point is computed using a plane fitted to its neighbors, according to the following procedure: Finding the k n nearest neighbor points: To accelerate the computational search of the k n nearest neighbor points a quadtree, i.e., a tree data structure in which each internal node has exactly four children, is used. For this, the implementation by Walker [15] is used, which is based on the work of Frisken and Perry [16].
Fitting of the plane: For every point w w where K i is a set of neighbor points indices, are defined. By computing the centroid c c c i of the neighbor points w w w j the 3 × 3 matrix A A A i can be calculated using Eq.
(2). A A A i is symmetric and positive semidefinite, thus the eigenvalue decomposition can be determined by Eq.
contains the corresponding eigenvectors.
Assuming that the eigenvalues are arranged in decreasing order, then n n n i ∶= ±v v v 3 i represents the normal vector of the fitted plane to the nearest neighbor points of w w w i . The sign of the vector is then chosen according to the machining orientation within the machine tool coordinate system: Following the definition in Fig. 2 the Z-component of the normal vector is chosen to be positive.
In this way, for the point cloud P a set of normal vectors N P , one for each point w w w i , is obtained according to Eq. (4).
It may be noticed, that this computation of the normal vectors for a point cloud is already implemented in some commercial software, e.g., in the Computer Vision System Relations between a nominal surface profile f, its normal vectors n n n , the tool nose radius r T and the tool path Toolbox of MATLAB or open3D for Python. Since the computation of the normals plays a crucial role in the computation of the trajectory, in particular for big tool radii, also a modified version of [14] and additionally the method from [17] are considered. The latter is based on a weighted average of the normals to the neighboring triangles obtained by a Delaunay triangulation. However, as it will become clear in Sect. 3, using an equally distributed point cloud leads to almost optimal accuracy 1 for a suitable point density.

Calculation of tool center surface
Once the normal vector of each surface point is calculated, the surface on which the tool center point moves can be computed. The most common cutting tools used in diamond turning, especially for machining freeforms, have a circular cutting edge defined by the nose radius r T and the geometrical center point coordinates c c c T = (x T , y T , z T ) in the machine tool coordinate system. The tool modeling assumes an ideally sharp cutting edge, which together with the tool center point defines the rake face plane E T with the corresponding normal vector n n n T according to the definition given by where u u u are the points on the plane and the vector n n n T in the surface reference system is given by Here and C are the rake and the C-Axis angle, respectively. For the general situation, the tool rake face is not always parallel to the XZ-plane but tilted by the tool's rake angle.
For a given point w w w i ∈ P , i.e., a point on the surface to be machined, the corresponding tool center point is computed by projecting the surface normal vector n n n i onto the plane E T , accordingly scaling its length and adding the resulting vector ñ n n i to the surface point. The according mathematical procedure is described in the following and corresponds to the depiction in Fig. 4.
The surface normal vector component in direction of n n n T is subtracted from n n n i and scaled according to Eq. (7).
(7) n n n i ∶=n n n i ‖n n n i ‖ ,n n n i ∶= n n n i − ⟨n n n i , n n n T ⟩ ⋅ n n n T .
Finally, the surface (or the point cloud) P T on which the tool center point moves is computed by extending the surface points w w w i along the obtained vectors according the definition given in Eq. (8).

Calculation of the trajectory
Usually, the basis tool trajectory in face turning of optics is supposed to be an Archimedean spiral defined in the XYplane. Equation (9) defines the spiral T 2D in the Cartesian machine tool coordinate system in function of the time t using the cylindrical coordinates consisting of a parametrization of the C-axis angle and the radial coordinate r. The radial coordinate r corresponds to the distance to the spindle axis C.
Once P T is determined, the 2D trajectory is projected onto it along the Z-direction. By interpolating the Z-values within P T the tool path points are generated. The position of the trajectory points in the XY-plane remains unchanged during this operation. As a byproduct, using the mapping v v v i ↦ w w w i the cutter contact points are obtained.
By choosing the type of interpolation it is possible to improve the accuracy, see e.g., Theorem 4.4.20 in [18], where estimations of the interpolation error depending on the differentiability of the surface, the degree of the interpolation polynomial and the size of the mesh are given. Since for every point the neighbor points were already determined in order to compute the normals (see Sect. 2.1), it is further possible to interpolate only locally to improve the memory management during the computations.

Analysis of accuracy loss
In this section, the different error sources that affect the tool path accuracy loss as defined in Fig. 1 are analyzed for various surfaces. In a first study, only the calculation of the normal vectors for a point cloud, see Sect. 2.1, is considered and the results are compared with the exact normal vectors. The influence of the definition of the point cloud P, the number of nearest neighbors k n and the used methods are discussed. In a second study, the interpolation error when interpolating the Z-values of the tool path points onto the tool center surface P T is analyzed and the trajectories resulting from a point cloud and from the corresponding analytical calculation are compared.

Accuracy of the normals calculation
First, a freeform surface defined by Eq. (10) is considered, see also Fig. 5.
The study investigates the accuracy of the surface normals calculation, is conducted using b = 3 and is independent of the machine tool and the cutting tool. It is assumed that the unit of the coordinates is mm. The normal vectors are calculated for points (x i , y i ) lying within a circle of radius 10 mm, for a number of points k ∈ {1E6, 4E6} and different numbers of neighbor points k n between 3 and 20, for both a Cartesian mesh and quasi-random points from the Halton sequence. The first grid distribution is motivated by the fact that some measuring devices, e.g., a Fizeau interferometer, return values on a Cartesian mesh, while the second distribution because of its low discrepancy, i.e., the points are random but evenly distributed. Then, motivated by on-machine metrology two additional methods are considered, both relying on surface data from a spiral distribution, as if it were generated by a spiral scanning measurement. Taking an exemplary but usual configuration, the spiral points have a maximal arc distance ds of 8.7 m or 10 • angular distance, while the distance dr between every C-axis turn is set to be 8.7 m, resulting approximatively 4E6 points lying within a circle of radius 10mm. Once the data is interpolated on a Cartesian grid the new grid points are used to compute the normal vectors. In the second case the method from [17] is directly applied to the spiral points to obtain the normal vectors. In this case, the number of neighbors is determined automatically for each point (the user cannot control it) using a Delaunay triangulation and each normal vector is computed using a weighted average of the normals to the faces surrounding the point.
To calculate the error of the used methods for computing the surface normals, the angle e ang i between the estimated normal vector n n n i and the analytical formula n n n(x i , y i ) is calculated according to Eq. (11).
The resulting mean values of the error angles fork ∈ {1E6, 4E6} in function of the number of neighbor points k n are shown in Fig. 6. The main observation is that in the case of a Cartesian mesh the mean error has local minima at k n = 5, 9 and 13, which show an improved error behavior in comparison to all other configurations. In the way the search for the neighbor points is defined, the actual point is considered as its own neighbor, meaning that the real numbers of neighbors for the above-mentioned amounts are 4, 8 and 12. The resulting neighbor points configurations within the regular grid are depicted in Fig. 7. For these special configurations, independently of the surface, the minima are a consequence of the approximation property of the central finite difference, which is of second order with respect to the mesh width, see Appendix 1.
The approximation of a function's derivative depends on the distance of the points used to compute it. Thus, by increasing the distance of the neighbor points, the accuracy of the proposed method may decrease, as it can be seen by comparing the curves with 1E6 and 4E6 points in Fig. 6. In contrast to the Cartesian grid, the quasi-random mesh (Halton set) does not show these local minima because the second-order approximation property is not valid for quasi-random points.
In case of the interpolated spiral data, the local minima are observed again, with the difference that the best approximation is for k n = 13 points. This level of accuracy is reached also by the method from [17] applied directly to the spiral points, as shown in Table 1. In this first attempt to study the approximation of normal vectors it is noted that the choice of the grid points play an essential role in the determination of the normal vectors. First, the best results are obtained on a Cartesian mesh with suitable value k n . Second, in the case of the spiral it turned out during the computations that the spiral points should be chosen to be well distributed, namely the arc length between the points should be approximatively equal to the distance of the spiral turns, as it can be observed in Table 1. The case of the table's last row is computed to show how the form of the triangles affects the precision of the method from [17].

Accuracy of the tool path points
Since not only the angle error resulting from the surface normals calculation contributes to the tool path accuracy loss, also the contribution of the interpolation of the tool trajectory points needs to be taken into account. Therefore, this section considers the error sum due to the interplay of approximation of the normal vectors and the interpolation in order to quantify the tool path accuracy loss according to Fig. 1.
The trajectory for this investigation is an Archimedean spiral with a maximal radius R out = 10 mm where the spiral points have a maximal arc distance of 1 mm or 1 • angular distance, while the distance between every C-axis turn is set to be 10 m. The cutting tool is supposed to have a nose radius of 1 mm and a rake angle of 0 • . The XY-coordinates of the point cloud lie on a regular Cartesian grid. To compute the surface normals 5 neighbor points according to Fig. 7 are used.
For this study two surfaces are considered. The first is defined by Eq. (12), see Fig. 8. The second is defined by z 1 according to Eq. (10).   Fig. 7 Configurations to calculate the surface normals by fitting a plane to various chosen numbers of neighbor points k n The first surface z 0 has zero curvature and the tool center trajectory is computable by analytic methods and therefore considered to be known. Since it is a plane, the surface normals calculation does not lead to any inaccuracies, allowing to quantify the error contribution due to interpolation only for several angles . This case needs particular attention, because since the surface z 0 is tilted by it cannot be machined in the center region without further attention. Therefore, the center region of r < 0.5 mm is not considered by this investigation. The second is a general freeform surface, but since the exact solution for the tool trajectory points using the steady-X method cannot be determined by using only analytic mathematics, it is proceeded in the following way: It is supposed that the tool center point moves on the surface according to Eq. (10) and a point cloud P T is computed using Eq. (10). Then the point cloud P, defining the surface to be machined, is computed using the tool center surface (given by P T ). In this way the tool center surface P T is analytically defined by the surface Eq. (10). Using the points P, the point cloud P T can then be calculated according to the presented procedure. To evaluate the error the tool trajectory points interpolated within P T are compared to the exact values computed using Eq. (10).
To quantify the error of the trajectory, first the spiral points (x , y ) , = 1, … , N are computed, then the interpolation of the points in P T is used to obtain the values of the Z-axis Z . These are then compared with the exact values z = z(x , y ) , where z is the surface to be machined and chosen according to the Eqs. (13) and (14).
-Planar surface z 0 -Freeform surface z 1 By computing the two quantities the differences in Z-direction between the trajectory computed by the proposed method and the exact trajectory are obtained.
In the Figs. 9 and 10 the local errors e i are plotted in function of the coordinate in C-direction over one rotation for different inclination angles and two different radial positions r. A sinusoidal behavior of the error with the amplitude gradually increasing with the tool moving closer to the rotation center can be observed. This is supposed to be related to the decreasing wavelength of the sinusoidal changes of the Z-values when machining the inclined planar surface towards the center. The global error values depending on the inclination angle are given in Table 2, where, as expected, increasing error values for a more inclined planar surface are observed. The given maximal values e max are throughout < 4 nm and except the maximal value of 3.9002 nm even < 1 nm. However, it should be noted that with an increasing inclination angle also other difficulties when machining the surface, for example a sufficient tool clearance angle, arise.
To study the approximation properties for the surface z 1 it is proceeded as already explained at the beginning of this section. Starting from the definition that z 1,tool represents the tool center surface, the surface to be machined is computed by using analytical formulas to calculate the surface normals. Based on that, the tool center surface can be computed once more using the numerical method presented in the Sects. 2.1 and 2.2. This way, two point clouds for the tool center surface are derived, one by using exact mathematics and the other by using the presented numerical method. In addition, it should be noted that the analyzed error values may be greater than usual due to the repeated computation ig. 9 Using 1E6 surface points for a surface size of R out = 10 mm : Error between analytical and computed trajectory using the proposed method in the case of the planar surface z 0 with inclination of 2 • , 5 • and 10 • evaluated for a single spiral loop at the radii 2 mm (a) and 0.5 mm (b). The numerical experiment shows the dependence of the error on the slope and on the radius from the spindle axis ig. 10 Using 4E6 surface points for a surface size of R out = 10 mm: Error between analytical and computed trajectory using the proposed method in the case of a planar surface z 0 with inclination of 2 • , 5 • , and 10 • evaluated for a single spiral loop at the radii 2 mm (a) and 0.5 mm (b). Again the dependence of the error on the slope and on the radius can be observed Table 2 Comparison of the error values for the surface z 0 (inclined planar surface) using a regular grid with 4E6 surface points and linear interpolation. Due to a singularity, which appears near the center as a consequence of the constant surface slope, the error has not been considered for r < 0.  Table 3 for different values for the amplitude b. As can be seen, the error values e max grow linearly as a function of the used number of surface points, which can be justified with the application of linear interpolation. Because of high slope angles and a certain variation of the Z-values depending on the amplitude b, a higher number of surface points is needed to achieve a reasonable small maximal error value along the trajectory. The surface with an amplitude of 3 mm ( b = 3 ) is viewed rather as an example to prove the accuracy of the proposed method, since the maximal surface slope angle exceeds the typical available cutting tool clearance angles. Even for the surface with an amplitude of 1 mm ( b = 1 ) the maximal error value can be reduced to about 1 nm by using a number of surface points of 16E6. All of the previous calculations are performed using linear interpolation. However, it is possible to achieve greater accuracy by increasing the interpolation order. To illustrate the changes in the error values due to a higher interpolation order, the computation and analysis for the surface z 1 with b = 3 are performed additionally applying cubic interpolation. The results given in Table 4 show a significant improvement by means of greatly reduced error values to almost 20 % compared to the values given in Table 3 using linear interpolation. It can be emphasized that even for the chosen case of b = 3 the result of e max is of sub-nanometer range for 16E6 surface points.
In Sect. 3.1 it is observed how much the method used to determine the normals affects the angle error defined by Eq. (11). To quantify their influence on the final trajectory the surface z 1 with b = 3 is considered with normal vectors obtained with different computation methods. To better show the effect of the normals cubic interpolation is used. As it can be seen in Table 5, the trajectory resulting from the quasi-random point cloud (Halton set) is at least two orders of magnitude less accurate than that computed with the Cartesian grid, while with the grid interpolated from a spiral distribution a similar maximal error is obtained. In the latter case, the mean error is increased remarkably, but remains below 1 nm. Also the tool path computed with the method from [17] shows an increased mean error, but still below 1 nm. The maximal error is slightly better than for the Cartesian grid. Investigations identified that this is caused by the double computation of the tool center points to obtain the exact solution, as explained at the beginning of this section. However, the standard deviation of the error as well as the mean error are smaller in the case of the Cartesian grid.

Analysis of machine tool interpolation
According to Fig. 1 the accuracy loss between the surface definition and the machined surface can be divided into two parts. The first is caused throughout the tool path generation and the second consists of the deviation resulting from the interpolation the machine tool controller applies on the programmed points as well as of the following error occurring in the position control loop of the axes drives. Whereas the first is analyzed and quantified in the previous Sects. 3.1 and 3.2, the second part addressing the accuracy loss related to the execution of the tool path points by the machine tool is not addressed so far. Hence, it is necessary to assess how the deviations of the tool path points compared to the exact calculation affect the execution by the machine tool. Although  the execution quality is supposed to be also influenced by the individual employed diamond turning machine tool as well as chosen parameters (NC point resolution, cutting speed, feed rates), the following analysis is assumed to represent common conditions of diamond turning. The analysis is divided into two parts, whereas the first investigates the tool path of the surface z 0 and the second the one of z 1 . Both utilize a diamond turning machine of the type LT-Ultra MTC 650 having a POWER PMAC control.

Analysis of the surface z 0
For the first investigation Z-axis positions following a sine function with the C-axis positions as arguments are programmed to the machine tool and executed. Subsequently, the resulting path interpolated by the control is analyzed. A nearly sinusoidal movement in Z-direction in function of the C-axis angle corresponds to machining an inclined planar surface (see geometry z 0 according to Eq. (12)). However, for the present study the small differences to a real sine due to the radial feed rate in case of a spiral machining path and the tool nose radius compensation are neglected. The tested sine wave has a period of 4s and an amplitude of 15 tan(5 • ) mm corresponding to a planar surface with an inclination of 5 • at a radius of 15 mm.
To investigate the effect caused by the numeric inaccuracies elaborated in Sect. 3.2, the Z-values of the programmed sine are overlaid with a Gaussian noise and the resulting axes movements are recorded by reading out the position signal of the machine tool's scales. As interpolation mode the spline mode of the POWER PMAC control is chosen, which interpolates programmed points by generating non-rational B-splines and is popular for executing curved paths represented by a multitude of points. It has the advantage of smooth axes movements, but the interpolated path does not exactly go through the programmed points. Hence, the interpolated sine wave commanded to the axes drives shows a slightly smaller amplitude compared to the initially programmed points. Figure 11 depicts one sine period consisting of points generated by the machine tool control based on programmed points without an overlaid noise and its deviations to a least squares sine fit. To exclude the effect of the reduced amplitude as far as possible the fitted sine function is supposed to represent the nominal movement for the analysis. This assumption is confirmed by the maximum deviation of just about 2E−4nm shown in Fig. 11.
The small standard deviations stated in Table 2 being in the sub-nm range are expected to result by the fact that the normals calculation for the flat surface z 0 is free of numerical inaccuracies. In this context it is expected that values below 1 nm barely affect the tool path execution. Therefore, the standard deviation n of the overlaid Gaussian noises are chosen according to the values for the freeform surface z 1 given in Tables 3 and 4. With the value of 3.5nm all values of Table 4 are covered, while the value of 5nm covers most of the values in Table 3. Their effect is determined by the difference of the resulting movement path and the nominal sine function. Figure 12 shows the resulting movement errors over one sine period, which Commanded Z-value Deviation to fitted sine function standard deviations lie within a range of about 1nm and are 5.8138nm (no overlayed noise; Fig. 12a), 6.4822nm ( n = 3.5nm ; Fig. 12b) and 6.7533nm ( n = 5nm ; Fig. 12c).
Considering these results, it can be observed that the overlaid Gaussian noises barely affect the movement of the machine tool's Z-axis, since the standard deviations of the movement errors increase only by about 1nm compared to no overlaid noise. Due to the fact that the computation accuracy of the presented method can be increased even to values of e max = e mean = < 1nm (see Tables 2, 3 and  4), the movement error and therefore also the tool path accuracy loss according to Fig. 1 resulting out of numeric inaccuracies can be assumed to be negligible.

Analysis of the surface z 1
The second investigation addresses the tool path of the freeform surface z 1 by analyzing the tool path execution of one C-axis rotation at a radius of as shown by the exact Z-values in Fig. 13. Therefore, the exact solution and the tool path calculated by the proposed method are executed on the diamond turning machine while the actual axes positions are recorded from the axes' scales. By calculating the differences of the Z-coordinates between the two path executions the influence of the numeric inaccuracies is evaluated. Similar to the analysis in Sect. 3.2 the maximum error e Mmax , the mean e Mmean and the corresponding standard deviation M over all recorded position points are calculated. The NC points programmed to the machine tool control have a spacing in C-direction of 0.1 • and a feed rate in X-direction per rotation of 10 m.
The procedure is carried out for all parameter sets included in Tables 3 (tool center points determined by linear interpolation) and 4 (tool center points determined by cubic interpolation). Additionally, each parameter set is carried out with two numbers of C-axis rotations per minute ( 20 rpm , 40 rpm ) to study the influence of different path execution velocities. Besides the exact tool path, Fig. 13 also shows the differences over one C-axis rotation for the parameter set of b = 3 , 16E6 points for the tool path generation using linear interpolation and a C-axis rotation speed of 20 rpm. Figure 14 shows the resulting path execution errors corresponding to the parameters listed in Table 3 for the case of linear tool path interpolation and Table 4 Tables 3 and 4. However, towards their minimum values in Fig. 14 they reach a relatively constant level, whereas the tool path errors (Tables 3 and 4) further decrease to even sub-nm values. It is supposed that this indicates the limit of the machine tool's motion accuracy, which is also confirmed by the standard deviations of the movement errors stated in Fig. 12 with about 6 nm being of the same order of magnitude as the minimum values in Fig. 14. The maximum error values e Mmax do not fall below 20nm, which is explained to be caused by the outlier data points being part of the recorded position signal from the machine tool's scale as shown in Figs. 12 and 13. Furthermore, the experiments conducted with 40 rpm C-axis rotation speed show slightly lower error values for most of the tested parameter sets. This can be interpreted in the way that the increased motion accelerations to some extent inhibit the axes to follow the inaccuracies of the tool path calculated using the proposed method. Thus, the executed motion shows a slightly better agreement with the execution of the exact tool path with 40 rpm than 20 rpm.

Corrective machining results
The direct advantage of the proposed method lies in the corrective machining using point cloud based measurement data. Since the tool path generation is also point cloud based, no surface fitting operation needs to be applied to the measured data and it can be easily compared to the nominal surface definition by calculating the difference in Z-direction. Prior to this, a 2D robust Gaussian regression filter is applied to remove noise and outliers.
To test this functionality a tilted flat according to Fig. 8, whose Z-values follow the function z F = x tan(5 • ) , is machined on a diameter of 20 mm in oxygen-free copper using a spiral tool path towards the center (X = 0). This surface geometry is chosen, because its tilt requires three-axes machining similar to freeform surfaces while the flat's form deviations can be verified using an external laser interferometer, which unlike a freeform metrology device was available for this work. Figure 15 depicts the machining setup for the experiment utilizing a diamond turning machine of the type LT-Ultra MTC 650, a monocrystalline diamond cutting tool with a nose radius of 0.6 mm and an interferometric 1D distance probe system. Table 6 contains the parameters used for the tool path generation, machining process and on-machine measuring process. The number of surface points 4E6 and the linear interpolation mode for the numeric tool path generation are chosen, because this setting results in small motion errors near the obtained minimum values for the surface z 1 with b = 1 while requiring less computation effort than with using 16E6 surface points (see Fig. 14). Since the surface z 1 with b = 1 has stronger curvatures than the tilted flat used for the corrective machining experiment, even smaller motion errors can be expected, which justifies the chosen parameters.
One corrective pass is carried out based on an on-machine surface form measurement using the optical distance probe. Figure 16 shows the form error results of the first (Fig. 16a; without correction) and second ( Fig. 16b; with correction) machining pass. The peak-to-valley form error has been reduced from 1.195 mm to 0.467 mm, which corresponds to a reduction of 61% . An analysis of the initial form error (Fig. 16a) results in a cutting tool offset in X-direction superimposed by a temperature deflection in Z-direction being the main reasons for the pattern of the form error. The correction   Fig. 15 Experimental setup for the corrective machining of the tilted flat of the characteristic step at the center (at X = 0; shown in the sectional view of Fig. 16a), which, as also stated by Preuss [19], is the result of a tool offset when machining a tilted flat, is impeccable. However, there seems to be a residual impact of a temperature deflection, which has not been part of the systematic errors influencing both machining passes in the same way. Besides, the surface's roughness after the corrective pass of Sa = 3.0nm (according to ISO 25178 and applying a Gaussian filter with cut-off wavelength of 80 m is of typical order of magnitude for diamond turning.
Carried out investigations indicate that the outlet temperature of the cooling mist nozzle causes the temperature deflections during the machining passes by cooling down the workpiece when streaming against it during machining. Figure 17 depicts recorded data showing that the mist outlet temperature is about 0.3 • C lower than the ambient temperature. The dashed line part shows the temperature measured 20 mm in front of the mist nozzle while the mist is turned off, whereas the solid line part while the mist is turned on like it is during machining. Figure 18 shows the relations between a decreasing temperature during the machining process from the circumference to the center and the resulting workpiece surface profile. For the sake of simplicity both profiles are depicted as linear function. The initial machining result is shown in Fig. 18a, where a temperature decrease ΔT 1 leads to a positive surface form error ΔZ W1a at the center, which is consistent with Fig. 16a. Supposing the same temperature profile ΔT 1 would occur during the corrective machining the surface form error profile would be compensated to the horizontal profile ΔZ W1b . Because of slow temperature assimilations a more moderate decrease ΔT 2 during the corrective machining can be assumed. This together with the correction based on the profile ΔZ W1a leads to the final surface profile ΔZ 2 , which corresponds to an over-correction.
The achieved form error reduction of 61% agrees with the reported values (between 57% and 61% ) recently published by Tong et al. [20], who focus on closed-loop freeform machining. The peak-to-valley values presented there lie between 207nm and 384nm and are achieved over half of the diameter considered in the present work using a fast tool servo, implying a reduced machining duration compared to the present work. This, together with the identified thermal deflections, indicates potential for further improvements. Thus, a more detailed analysis of the process and the employed on-machine metrology setup including its calibration and an uncertainty analysis will be focused on in a subsequent publication.

Conclusions and outlook
A new method for the computation of the cutting tool path in ultra-precision diamond turning, especially for threedimensional freeform and corrective machining using form measurement data, is presented and analyzed. It allows the handling of surfaces defined by equations and by point clouds without using symbolic computation methods. Instead, numeric tools are used for the computation of the surface normal vectors as well as the surface on which the center point of the cutting tool's nose radius moves. Together with interpolation techniques, the generation of a spiral tool path in case of the steady-X method ensuring a uniform axis motion in X-direction is provided. By a comprehensive analysis, the errors occurring throughout the computation steps are compared to the exact solutions for various geometries with diverse curvatures. Under realistic assumptions regarding the curvature of the surface to be machined and reasonable refinement of the applied numeric methods the results show errors < 1nm . Combined with an investigation of the effect of non-exact tool path points on the behavior of an exemplary diamond turning machine this allows to neglect the accuracy loss compared to the exact tool path. Besides providing a tool path of sufficient accuracy, the method offers to update the point cloud representing the surface to be machined by surface form measurement data. Presented form error results of corrective machining a flat, which is tilted by 5 • with regard to the XY-plane, showing a reduction of the peak-to-valley form error of 61% verify the ability of the proposed method for this application.
The results show the subtle interplay of grid, normal vectors and interpolation in the accuracy of the computed trajectory. While the role of interpolation is well known, additional points need to be addressed, e.g., the point distribution and the influence of noise in the computation of the normals. To the best of the authors' knowledge, even in recent works, where point clouds are used, this is not investigated.

Appendix 1: Analysis of the plane fitting method
From the results depicted in Fig. 6, it is observed that for point symmetric stencils on a Cartesian grid, the angle error behaves a lot better than for other general stencils. This suggests that it is depending on the surface, but this is because the plane approximation for point symmetric stencils delivers a second-order approximation of the normal vector. The aim of this appendix is to prove this claim for a smooth surface, or at least a three times continuously differentiable surface. Without loss of generality, the following configuration for a small number h > 0 is considered under the assumption that Then it follows where the missing entries are completed by symmetry. By letting A A A = X X X ⋅ X X X T = 2h 2