Lines of Curvature for Log Aesthetic Surfaces Characteristics Investigation

: Lines of curvatures (LoCs) are curves on a surface that are derived from the ﬁrst and second fundamental forms, and have been used for shaping various types of surface. In this paper, we investigated the LoCs of two types of log aesthetic (LA) surfaces; i.e., LA surfaces of revolution and LA swept surfaces. These surfaces are generated with log aesthetic curves (LAC) which comprise various families of curves governed by α . First, since it is impossible to derive the LoCs analytically, we have implemented the LoC computation numerically using the Central Processing Unit (CPU) and General Processing Unit (GPU). The results showed a signiﬁcant speed up with the latter. Next, we investigated the curvature distributions of the derived LoCs using a Logarithmic Curvature Graph (LCG). In conclusion, the LoCs of LA surface of revolutions are indeed the duplicates of their original proﬁle curves. However, the LoCs of LA swept surfaces are LACs of different shapes. The exception to this is when this type of surface possesses LoCs in the form of circle involutes.


Introduction
Curve and surface design for shipbuilding started during Roman times [1]. The progress of representing this in the form of mathematical equations and algorithms increased drastically with the introduction of numerical control (NC) to handle milling machines in the 1950s. In the 1960s, automobile and aircraft designs were the two key fields employing mathematicians to work on curves used for effective and efficient design.
Two classical curves used for computer aided design (CAD) and computer aided manufacturing (CAM) are Bezier curves and B-Splines [2]. Later, Nonuniform Rational B-Spline (NURBS) became the de facto standard for CAD/CAM since it can represent conics analytically [3]. However, these curves cannot be used directly for manufacturing due to their complex forms of curvature; hence, various fairing algorithms were devised to reduce the oscillation in their curvature profile [4][5][6]. Spirals that have monotonic curvature profiles became a good candidate for use in aesthetic design [7][8][9]. In 2005, Log Aesthetic Curve (LAC) was first introduced. This curve consists of spiral families, e.g., clothoid, logarithmic spiral, circle involute and Nielsen's spiral [10].
LAC was first proposed by Miura [10], by representing a Logarithmic Curvature of Graph (LCG) as a linear equation with α as its gradient. While the curvature profile involves the second derivates of curves, LCG involves the third derivates, hence making it suitable as a higher order shape interrogation tool [11]. The LAC family has a monotonic curvature profile which can be used to represent well known spirals, e.g., clothoid, logarithmic spiral, circle involute and Nielsen's spiral. Research in LAC started to progress rapidly after Yoshida and Saito [12] proposed a G 1 algorithm to interactively render LACs that are

Log Aesthetic Curves
A linear representation of a Logarithmic Curvature Graph (LCG) is the fundamental equation to describe LACs. The equation for LCG which has a slope, α, is shown below [10]: where s is the arc-length of LAC, ρ L is radius of curvature and C is a constant. After differentiation and simplification of Equation (1), we obtain: where Λ is the extra degree of freedom of LAC and 0 < Λ < ∞. Upon integrating Equation (2), we have: where ρ 0 is the initial radius of curvature of LAC, and ρ L varies from 0 to ∞. The arc-length, s, has an upper boundary and lower boundary which depends on the α. The curvature of LAC is shown below: By setting torsion, τ(s) = 0 in the case of planar curve, we can compute LACs by solving the Frenet-Serret formula: Table 1 shows the upper boundary and lower boundary of s with respect to α derived by Yoshida & Saito [12]. Since ρ 0 α ∈ [0, ∞], we must choose a suitable value of α in order to obtain a well-defined LAC segment.
To generate LoCs on the surface, we need to solve the IVP below numerically: where κ p is the principal curvature representing either the maximum principal or the minimum principal curvature. The non-zero factor (η and µ) can be obtained by normalizing the first fundamental equation as shown below: The sign of dt ds can be determined from Equation (14) [22]. LoCs can be generated by solving the IVP (Equation (12) and (13)) using numerical methods such as the Runge-Kutta method and substituting their values into the original surface function. To investigate the characteristics of LoCs, Joo et al. [22] proposed a novel method to compute the curvature and torsion of LoCs. In order to calculate the LoCs, we need first to compute the geodesic curvature of LoCs. Equation (15) represents the proposed method to compute u , v and the geodesic curvature, κ g .
where U = N × t = N × c (s), The curvature can be calculated using the principal curvature and geodesic curvature, as shown below: Furthermore, the torsion and first order curvature's derivative of curvature on the surface can be computed using the following equation: where t = c (s),

CUDA Programming in Wolfram Mathematica
In this paper, we used parallel programming to speed up the computation time for generating LoCs on the surface. The Mathematica code is a simplified procedure for CUDA implementation, as it internalizes the memory allocation, copies, processes and synchronizes. There are four steps towards CUDA implementation [27]: 1.
the CUDA code is cached to a binary file optimized for the GPU selection; 3.
the compiled code is cached to avoid future compilation; 4.
a check is performed as to whether the kernel exists in the compiled code.
Hence, the only thing the user needs to do is write the CUDA kernels as shown in the Appendix A. To note, all the codes must be placed into a string (CUDAKernelCode) so that it can be executed in Mathematica. Function func2( . . . ) is a device function which can be called from __device__ or __global__ function but cannot be called from Mathematica's language. Function func (float * sum, mint listSize) is a global function which can be called from and returned to CPU. The sum is returned in the form of array and mint is the Wolfram Language integer that is similar to int. The index is calculated based on the current thread to be used for execution. If (index < listSize) is to ensure that the output buffer is not overwritten. Function func2(index) is evaluated at index and stored at sum[index]. Next, CUDAFunctionLoad is used to load the program which passes the kernel code as a first argument, the kernel name to be run ("func") as a second argument, the function parameters as a third argument and the last argument is the block size. The result is stored in Kresult. CUDAMemoryAllocate is used to allocate CUDA's memory to store the result. The example shows that 1024 buffers are created to store the result. KResult[bufferK, 1024] shows how the kernel is executed and we can retrieve the result by using CUDAMemoryGet [bufferK]. The buffer must be freed up by the user manually by using CUDAMemoryUnload [bufferK]. The complete implementation of CUDA coding in Mathematica for this paper is readily available in GitHub [32].

Log Aesthetic Surface of Revolution
For profile curve C = { f (u), g(u)}, that is a parametric curve defined in xz-plane, the general equation of the surface of revolution, S Rs revolving around z-axis is: Assume that the radius from the center of the surface of revolution is m, then we have the following equation for surface of revolution as: For example, if C is a circle, then S T is a torus with radius m. The first and second fundamental equation's coefficients of this surface are: x(u) ..
y(u) = dy du , and .. y(u) = d 2 y du 2 , respectively. We can choose any profile curve to generate this surface. The Super-spiraloid produced in [4] also has these settings, if one were to choose super-spiral as their profile curve. any profile curve to generate this surface. The Super-spiraloid produced in [4] also has these settings, if one were to choose super-spiral as their profile curve. Figure 1 shows LA surfaces of revolution by setting LAC as its profile curve with various types of . The thick black line is the LAC profile curve, , used to generate these LA surfaces of revolution. Note that we set = 3, = 0,2 and = 0, /3 for the entire computation. Since = = 0, the LoCs on LA surface of revolution are meridians (red) and parallel curves (blue) on the surface and they are orthogonal to each other. To reiterate, the red LoCs are duplicates of the original LAC used as the profile curve and the blue LoCs are circles with different radius of curvature.

Characteristics of LoCs for LA Surface of Revolution
The numerical methods used to compute the LoCs in CPU are the classical Runge-Kutta method (RK4) [14] with fixed step-size 0.001 and the Dormand-Prince method (DP) [15]. For the GPU, we used the classical Runge-Kutta method (RK4) [14] with a fixed stepsize of 0.001. DP is faster for LAC computation when compared to RK4 [15], Gaussian Kronrod [14], and Incomplete Gamma Function [13]. However, the DP's implementation in GPU needs extra effort. Since RK4 is computationally more efficient, we used RK4 for the GPU comparison in this work. Table 2 compares the CPU and GPU performance metrics by calculating the computation time to generate 22 LoCs (11 red and 11 blue lines) on

Characteristics of LoCs for LA Surface of Revolution
The numerical methods used to compute the LoCs in CPU are the classical Runge-Kutta method (RK4) [14] with fixed step-size 0.001 and the Dormand-Prince method (DP) [15]. For the GPU, we used the classical Runge-Kutta method (RK4) [14] with a fixed step-size of 0.001. DP is faster for LAC computation when compared to RK4 [15], Gaussian Kronrod [14], and Incomplete Gamma Function [13]. However, the DP's implementation in GPU needs extra effort. Since RK4 is computationally more efficient, we used RK4 for the GPU comparison in this work. Table 2 compares the CPU and GPU performance metrics by calculating the computation time to generate 22 LoCs (11 red and 11 blue lines) on the surface. We repeated these computations five times and the average computation time in CPU and GPU are 167.257 (DP), 1188.91(RK) and 8.28301 (RK) seconds, respectively. It is apparent that GPU greatly reduces the computation time for LoC generation.  Figure 2 shows the characteristics of LoCs of LA surfaces of revolution that correspond to the various types of α in Table 2 Figure 2 shows the characteristics of LoCs of LA surfaces of revolution that correspond to the various types of in Table 2.
(a) (b) The red LoCs on the surface are identical to each other since these are the duplicates of the original LAC. Hence, they have the same curvature and torsion profile. Figure 2a,b show that the curvature of the LoCs is indeed the curvature of its respective LACs in xzplane. When = −1, the curvature of the LoCs is linear, hence its derivative is a constant, which has the same curvature distribution as a clothoid. Figure 3 depicts the LoCs on LA surface of revolution with = 2 as a square pipe to keep track of the torsion. Figure 3a indicates the Frenet-Serret frames at each point. Figure 3b shows a close-up view of the surface and Figure 3c represents the LoCs on the LA surface of revolution. Based on Figure 3b, there is no twist along the LoCs, which is consistent with the fact that each curvature line has 0 torsion. The red LoCs on the surface are identical to each other since these are the duplicates of the original LAC. Hence, they have the same curvature and torsion profile. Figure 2a,b show that the curvature of the LoCs is indeed the curvature of its respective LACs in xz-plane. When α = −1, the curvature of the LoCs is linear, hence its derivative is a constant, which has the same curvature distribution as a clothoid. Figure 3 depicts the LoCs on LA surface of revolution with α = 2 as a square pipe to keep track of the torsion. Figure 3a indicates the Frenet-Serret frames at each point. Figure 3b shows a close-up view of the surface and Figure 3c represents the LoCs on the LA surface of revolution. Based on Figure 3b, there is no twist along the LoCs, which is consistent with the fact that each curvature line has 0 torsion.  Figure 2 shows the characteristics of LoCs of LA surfaces of revolution that correspond to the various types of in Table 2.

CPU(DP) CPU(RK4) GPU(RK4)
(a) (b) The red LoCs on the surface are identical to each other since these are the duplicates of the original LAC. Hence, they have the same curvature and torsion profile. Figure 2a,b show that the curvature of the LoCs is indeed the curvature of its respective LACs in xzplane. When = −1, the curvature of the LoCs is linear, hence its derivative is a constant, which has the same curvature distribution as a clothoid. Figure 3 depicts the LoCs on LA surface of revolution with = 2 as a square pipe to keep track of the torsion. Figure 3a indicates the Frenet-Serret frames at each point. Figure 3b shows a close-up view of the surface and Figure 3c represents the LoCs on the LA surface of revolution. Based on Figure 3b, there is no twist along the LoCs, which is consistent with the fact that each curvature line has 0 torsion.

Log Aesthetic Swept Surface
In this section, the characteristics of LoCs of LA swept surface are investigated. First, a LAC is swept along another LAC to form a complete LA swept surface. The general equation of the swept surface is obtained as we sweep along the path curve, ( ), with the trajectory of the profile curve, ( ) [33]:

Log Aesthetic Swept Surface
In this section, the characteristics of LoCs of LA swept surface are investigated. First, a LAC is swept along another LAC to form a complete LA swept surface. The general equation of the swept surface is obtained as we sweep along the path curve, C 1 (u), with the trajectory of the profile curve, C 2 (v) [33]: where T 1 is a rotation matrix. Suppose that the path curve, C 1 (u) = {x 1 (u), y 1 (u), 0} and profile curve, C 2 (v) = {0, −x 2 (v), y 2 (v)} are LACs, then we sweep curve C 2 (v) along curve C 1 (u) to form a complete LA surface. The equation of the LA surface can be computed by applying the details into Equation (20) to form: Since . C 1 (u) is a normalized tangent, then we have t 11 (u) 2 + t 12 (u) 2 = 1. Hence, the LA surface can be simplified as follows: Next, the first and second fundamental equation's coefficients of LA swept surface are: . t 12 (u) = 0, .. ..
ized tangent of C 2 (v). Take note that F = M = 0, indicating that the LoCs are similar to parametric lines.
As a numerical example, a LA surface is generated by setting α 1 = 2, and Λ 1 = 0.1. As shown in Figure 4, the initial radius of curvature, ρ 10 = 1 for curve C 1 (u) and α 2 = 2, followed by Λ 2 = 0.1 and initial radius of curvature, ρ 20 = 0.5 for curve C 2 (v). Based on Figure 4b, the green curve represents curve C 1 (u) while the black curve represents curve C 2 (v). As for Figure 4c, the zebra map of the LA surface indicates a smooth surface. parametric lines.
As a numerical example, a LA surface is generated by setting = 2, and Λ = 0.1. Figure 4, the initial radius of curvature, = 1 for curve ( ) and = 2, followed by Λ = 0.1 and initial radius of curvature, = 0.5 for curve ( ) .Based on Figure 4,b, the green curve represents curve ( ) while the black curve represents curve ( ). As for Figure 4c, the zebra map of the LA surface indicates a smooth surface.

Characteristics of LoCs for LA Swept Surface
In this section, we investigate one part of the LA swept surface as shown in Figure  4d. The red lines (labeled as to ) are the LoCs in maximum principal direction and the blue lines (labeled as to ) are the LoCs in minimum principal direction. First, we compared the computation time of 12 LoCs (consisting of six lines in the maximum principle direction and six lines in the minimum principle direction) on the CPU and GPU. Figure 5 shows swept LA surfaces with different and its respective LoCs.

Characteristics of LoCs for LA Swept Surface
In this section, we investigate one part of the LA swept surface as shown in Figure 4d. The red lines (labeled as L 0 to L 5 ) are the LoCs in maximum principal direction and the blue lines (labeled as D 0 to D 5 ) are the LoCs in minimum principal direction. First, we compared the computation time of 12 LoCs (consisting of six lines in the maximum principle direction and six lines in the minimum principle direction) on the CPU and GPU. Figure 5 shows swept LA surfaces with different α and its respective LoCs.      Figure 7 shows the curvature (κ) and derivative of curvature (κ ) for LoCs in the maximum principal direction (L 0 to L 5 ) for LA swept surfaces, as illustrated in Figure 4d. It is apparent that the curvature profiles and their derivatives of LoCs are monotonic.
We further plotted LCG to verify that these LoCs are indeed LACs as shown in Figure 8. The LCG plots share similar colors with the LoC, as seen in Figure 7. The LCGs are all in the linear form with a fixed LCG gradient ∆ i as stated in the legend. Hence, all L 0 to L 5 LoCs on each LA swept surface are confirmed to be LACs. However, they are dissimilar to their LAC's original profile, except for the case of α = 2 which represents the family of circle involutes. For the rest of α, the LCG gradients vary gradually as we sweep along another LAC profile curve. Figure 9 depicts the curvature and its first derivative of the blue LoCs of LA swept surfaces as shown in Figure 4d. The LoCs (D 0 to D 5 ) of each LAC swept surface have the same curvature distribution as its original LAC profile curve. Figure 10 shows its LCG plot which verifies that each set of minimum principal direction of LoCs is indeed similar to its original LAC profile curve.

Conclusions and Future Work
We presented two types of LA surface along with their respective LoCs in this paper. In order to speed up the computation time, we implemented LoC computation on the GPU, which resulted in the GPU performing approximately 20 times faster as compared

Conclusions and Future Work
We presented two types of LA surface along with their respective LoCs in this paper. In order to speed up the computation time, we implemented LoC computation on the GPU, which resulted in the GPU performing approximately 20 times faster as compared to the CPU for LA surfaces of revolution. For LA swept surface, the GPU computation is 197 times faster than CPU. Next, we analysed the curvature and its derivatives of LoCs for both types of surfaces. As expected, the LoCs on a LA surface of revolution are also parametric curves in which the meridians are LACs of their profile curve, and the parallels are circular arcs. However, for LA swept surface, the curvature and its derivative of LoCs in the minimum directions are indeed LACs with original LAC profile. However, the LoCs for maximum principal directions are dissimilar except in the case of α = 2. This indicates that the LoCs and parametric curves of LAC swept surface for this case are all circle involutes that did not change due to the sweep. Therefore, in the case where the designer wants uniform LoCs of the LAC swept surface, this can be achieved by setting α = 2.
Our future work includes the implementation of LA space curves to generate a free form surface design, for example, using Coon's patch, and investigating its LoCs. We are also planning to extend this work by investigating the LoCs for surfaces developed with super-spirals [17].