Viewpoint Generation Using Geodesics and Associated Semi-Automated Coverage Path Planning of Panels for Inspection

Featured Application

The specific application of this work is in the robotic path planning of camera-based non-destructive testing systems such as active thermography.

Abstract

The coverage of a surface using multiple viewpoints is a topic of great interest for robotic path planning in inspection applications. Two approaches for coverage path planning are broadly addressed in the literature—geometric methods and optimization methods. While the optimization methods may be the most flexible, they are frequently difficult to implement in practical applications due to their NP hard nature. We present here a geometric algorithm for the coverage path planning of panels used for aerospace applications using a generic camera model that can represent area inspection techniques like thermography and laser shearography. This algorithm relies on drawing a 2D grid on the 3D surface of the panel using geodesic lines on the surface. The coverage of the surface is performed by propagating geodesic lines from a starting point until the patch thus covered diverges too much from a flat surface, and after that, the coverage is continued from another point. The propagation of the geodesic lines is stopped when they begin to converge or diverge, and we define two criteria for the stoppage. We show that the proposed algorithm has good results for 3D virtual models and emphasize its speed, simplicity, and reliability for such applications.

1. Introduction

In certain applications in robotics and autonomous vehicles, we desire to cover an area, and we would like to obtain a path by which the entire area is covered in the best way using some metric such as the shortest distance, lowest travel time, or energy expenditure. Examples of such applications include vacuum cleaner robots covering the floor, lawnmower robots cutting grass on a lawn, plowing robots covering a field, inspection of underwater oil rigs and pipelines using autonomous submarines, and so on. This generic problem is called coverage path planning. Ref. [1] gives an overview of the generic applications of coverage path planning.

In the case of inspection applications where the sensor is held at a distance, what we desire to know is a set of viewpoints from which the entire surface can be inspected using the sensor at those viewing points. This problem is generically known as the View Planning Problem (VPP) and in case we need to plan a route from one viewpoint to another (such that all viewpoints are covered), then that problem is a subset of the Traveling Salesman Problem (TSP). Both VPP and TSP are well known to be NP hard problems and their combination TVPP (Traveling View Planning Problem) is thus also NP hard [2]. There are also other ways of framing the problem such as the art gallery problem, the watchman route problem, and so on (see [1]) but the common solution to these is based on numerical optimization. Ref. [3] is another survey on coverage path planning algorithms, but this too covers mostly optimization-based techniques.

The problem addressed here is the imaging of a curved panel for the purpose of Non-Destructive Testing (NDT), so as to cover the entire surface. We assume that the imaging will be done by a device like a camera. The device (camera-based NDT system) will be moved to a set of viewpoints, oriented at the appropriate angle and then the images would be taken. Since the viewpoint can be placed anywhere in space and oriented in any direction, in general, there may be an infinite number of viewpoints from which the surface may be viewed. In our paper [4], we show a method for planning the coverage by manually selecting each viewpoint. Figure 1a shows a generic viewpoint at which a camera is placed and oriented with its viewing frustum capturing the image of the panel from that viewpoint (according to the parameters of the view camera). Considering that a surface mesh is a discrete representation of the panel surface consisting of surface elements like polygons, we can reduce the solution space by just assuming a single viewpoint for each surface element which is normal to the surface at that point. This is shown in Figure 1b, where we have generated a viewpoint for each surface patch (the surface patch here being one triangle of a triangle mesh) using the normal to the surface. Each viewpoint can capture more than one surface patch (in this example, viewpoint V6 covers surface patches S4 to S8). Then the optimization problem becomes a matter of choosing a minimum number of viewpoints such that all the surface patches are viewed.

If we use the above approach with the surface modeled as a triangle mesh we will be having thousands of surface elements and viewpoints even for a simple surface and these will all be variables for optimization and thus would need a commercial solver to solve it in any reasonable amount of time. Ref. [2] shows a typical approach to solving TVPP with numerical optimization using IBM Cplex solver, which is not only expensive but difficult to integrate into custom applications. Another example [5] explores the topic of coverage path planning of outdoor structures with UAVs using optimization-based methods, the results of which can be seen in Figure 2A. In [5], the church model has been simplified to just 200 triangles in order to keep the optimization algorithm from taking a massive amount of time. Ref. [6] explores coverage path planning for structural inspections using a drone, and various heuristic methods to solve the optimization problem (e.g., genetic algorithms ) are implemented here.

An alternative approach to these optimization methods and solving NP-hard problems is to use geometric algorithms; however, care must be taken to use the appropriate algorithm for the problem at hand. The main contribution of this paper is to demonstrate the usage of a geometric method instead of complex numerical optimization methods and show that we get good results for the 3D virtual models of real-world examples. Ref. [7] describes a geometric algorithm using slicing planes for ‘circling’ around an underwater structure using an underwater robot for inspection (see Figure 3a), and Ref. [1] also describes a similar circling strategy for UAVs in urban environment. Ref. [8] describes a basic geometric algorithm for coverage path planning for spraying using drones.

Figure 2. Coverage path planning of structures using UAVs. (A) UAV inspection of church using optimization methods [5]; (B) UAV inspection of wind turbine using slicing-based geometric method [9].

Figure 2. Coverage path planning of structures using UAVs. (A) UAV inspection of church using optimization methods [5]; (B) UAV inspection of wind turbine using slicing-based geometric method [9].

Geometric algorithms based on slicing the surface based on planes are commonly used in 3D printing applications; however, if the distance between planes is kept constant then we will immediately run into problems for even the simplest curved surfaces (see Figure 3b for how a gap in between viewpoints has occurred), and in general such a technique cannot be used efficiently without making some assumptions about the surface mesh. Ref. [7] have implemented the circling algorithm for their underwater inspection and note that it is ‘inflexible’ and ‘not optimal’. Slicing algorithms have also been implemented in [9] for the purpose of inspecting large structures like wind turbines using UAVs (Figure 2B). Although geometric algorithms do not guarantee optimality or complete coverage, it is possible to use the geometric methods for applications such as panel inspection for aerospace systems, and obtain efficient algorithms that give near-optimal solutions by selecting the appropriate input parameters carefully. In addition, geometric algorithms are well known for their speed, simplicity, and ease of use.

We proposed a simple geometric approach based on using the geodesics of the surface, in [10]. In this approach, the surface is first partitioned into near-developable patches. On each patch, the near rectangular camera image sensor projections are placed, side by side, with overlaps to ensure that all points on the surface are covered. What is done is to place one rectangle, and then construct adjacent rectangles by moving in the length and width directions along the geodesics of the surface in those directions. The advantages of using geodesics are that (a) on near planar surfaces geodesics can form a near rectangular grid, and (b) geodesics are well-defined curves that can be constructed using simple algorithms.

For finding geodesics, we use the approaches proposed by [11,12]. Ref. [11] gives an algorithm to compute a discrete geodesic path on a surface mesh, and Ref. [12] uses this algorithm to fit woven cloth (resembling a 2D grid) onto the surface. We proceed in a similar fashion as [12] but dispense with the notions of warp, weft, and other terminologies related to fabric draping and replace the node mapping algorithms in [12] with our own node insertion algorithm. We show that the algorithm achieves good results for real-world scenarios.

In [10], we have described the underlying process as follows. Figure 4a shows a generic shape that can be decomposed into several sections, each of which is near developable, i.e., a 2D grid can be mapped on the surface without much distortion. Going from the generic to the specific, this model of using sections as parts of surfaces is exactly what is used in the aerospace industry (see Figure 4b), and the surface of the aircraft or launch vehicle is made of a large number of such sections which are also called panels. These panels are usually manufactured (and inspected) separately and then fitted together into the final vehicle. We show our results on panels taken from the real-world aerospace industry—an octant that is part of a propellant tank and a cylindrical surface with cutouts that is part of a launch vehicle. The notion of constructing a non-developable shape out of sections is far from being restricted to the aerospace industry; see, for example, the construction of a volleyball in Figure 4c. Thus, our technique can be used for a number of other applications as well.

In [10] we did not address the problem of partitioning complex doubly curved surfaces into developable patches. No method was given as to how to split up a surface in a way so that each section is near developable, although some manual methods are explored (e.g., Figure 5b where the hemisphere is split into five sections and coverage path planning is performed on each section separately). The entire surface cannot reasonably be mapped by a single grid, or we will run into problems (such as the hemisphere in Figure 5a), because of the non-developable nature of the shape in general, but a section of the surface will always be near-developable.

In this paper we improve on our previous concept by introducing the notion of near-developability, i.e., at what point the node propagation algorithm must be stopped in order to avoid the converging or diverging of geodesic lines on the surface. Two approaches are presented here—the ‘Indirect’ approach, which uses the concept of the Gauss map of the surface normals, and the ’Direct’ approach which takes into account the convergence or divergence of a viewpoint node from its neighbors. The results are explored for a wide range of panel shapes relevant to aerospace applications. The coverage algorithm here is “Semi-automated”, in which the user is expected to select the starting criteria for node propagation and then save the viewpoints. This provides far more versatility than the earlier ’Fully automated’ mode which has several failure conditions.

The method presented in this paper has the advantage over optimization-based techniques in being fast and simple to implement. Its advantage over other geometry-based techniques is that it is much more versatile than slicing-based methods. Finally, it has an advantage over just manually planning the path (e.g., as presented in our paper [4]) in that it is semi-automatic and requires much less human intervention than planning the pose of each individual viewpoint.

It is to be noted that the current paper focuses primarily on planning the viewpoints given start criteria provided by the human user, by decomposing the surface into regularly arranged neighboring camera viewpoints. The user input is very important to achieving proper coverage using this method, thus making it a semi-automated technique. Once the set of these viewpoints has been generated, they can be linked together into a path. However, the travel planning process can be done trivially, or by using known methods to solve the traveling salesman problem (see Section 6). We also note that in the given application of camera-based non-destructive testing, the position and orientation of the viewpoints and the coverage of the surface that they provide is of paramount importance, while optimally planning the travel plan or sequence of viewpoints is not that important because the robot can move very fast between viewpoints, but needs to hold at a viewpoint for much longer time to record the data. Hence, the algorithms described in this paper are basically about creating a set of viewpoints for covering the surface.

The rest of the paper is organized as follows. Section 2 gives an overview of the user workflow and inputs needed for the entire process. Section 3 details about the discrete geodesic path algorithm that is used to generate curved paths on the surface, and Section 4 explains our node insertion algorithm that maps nodes all over the surface to generate viewpoints and also incorporates the criteria for which the node propagation needs to be stopped. Section 5 shows the results of the coverage algorithms on a variety of panel shapes. Section 6 discusses how to link together all the viewpoints into a motion path, and finally Section 7 concludes and summarizes the paper.

2. Overview of Process Workflow

In this section, we detail the overall workflow of the process of generating the coverage path on the panel. First, we explain how the surface model is loaded into the 3D viewer that has been developed in C++ and the OPENGL framework. Since this is a semi-automatic technique, the user will be prompted to enter the initial criteria for starting the node propagation process. The node propagation will continue as long as the stoppage criteria are not met. After the node propagation process is stopped, the user can save these viewpoints and start again at another point. They can continue this process until the coverage is completed.

2.1. Loading Model into 3D Viewer

The virtual models used in this effort can not be easily represented using MATLAB R2023b or any other easily available scientific simulation software. Instead, a customized 3D viewer was developed using OPENGL and C++. This software plays an integral role in all of the algorithms to be discussed in further sections. Figure 6a shows the virtual model of our panel as it is loaded into our 3D viewer. Figure 6b shows the actual panel being inspected (it is an octant shape that is part of a propellant tank, and the cutout is not included in the virtual model for simplicity). The triangle mesh (STL file) of the panel is loaded and placed in the 3D viewer. The blue region with white grid lines represents the floor, with the grid representing an area of approximately 2 m × 2 m. The virtual model of the panel is shown placed on the floor. The X, Y, and Z axes are represented using the red, green, and blue lines, respectively.

2.2. Selecting the Starting Point on the Surface

The starting point is selected by right-clicking on the surface. This point is translated from the 2D output frame (viewport) into 3D space and a ray is generated using ray casting techniques. The intersection of the ray with the surface is found by checking ray-triangle intersections with all the triangles in the mesh. If more than one intersection is found, we select the closest point. The algorithms for ray-casting and ray-triangle intersection are done using standard techniques common in computer graphics (refer, e.g., https://www.scratchapixel.com/, accessed on 26 December 2023 ) and are not repeated here.

The primary axes for the coverage algorithm (i.e., the X and Y directions along which the geodesic needs to be propagated) then need to be defined. By default, the primary axis of the 3D viewer (i.e., X-axis) is projected onto the tangent plane located at this point, and this defines the X or tangent direction. The Y direction or bi-normal vector can be found using the normal and tangent vectors. This process is shown in Figure 7a. The X and Y (or tangent and bi-normal directions) are used as the starting direction to project geodesics on the surface up to a certain distance. These geodesics are colored red (X direction) and green (Y direction). Another two sets of geodesics (colored yellow) are shown for diagonal directions. The algorithm for generating the geodesic paths is explained in more detail in Section 3.

2.3. Rotating the Starting Direction on the Surface

The user can rotate the start direction (i.e., X-axis or tangent vector) on the tangent plane of the starting point using the scroll button while holding down the right mouse button. An example of rotation of the start direction is shown in Figure 7b.

The start direction can have a significant impact on the node propagation on the surface. This can be seen in Figure 8b. We can select this in order to cover some types of surfaces more effectively.

2.4. Initiating Node Propagation

The user can send the command to initiate the node propagation process (the algorithm for this is described in more detail in Section 4). This creates a preview set of viewpoints using the nodes generated by this process. The camera model at the location corresponding to a viewpoint is shown in green (for current set of viewpoints), and yellow (for saved viewpoints). Note: in some figures in Section 4, it is shown in red. The outlines corresponding to these viewpoints are displayed on the panel. The outlines are colored cyan (for present viewpoints), and yellow-gray (for saved viewpoints). If the user is unhappy with the preview, they can select another starting point and starting direction and start again. Figure 8a shows the output after the preview is done.

2.5. Saving Viewpoints

If the user is satisfied with the preview (Figure 8a), they can save the current set of viewpoints by sending the appropriate command. Upon doing so, these viewpoints are saved and they appear in a different color (i.e., yellow) on the panel (the outlines and camera locations). After this, the user can select another point (Figure 9a and start the preview from there as well (see Figure 9b). They can save those viewpoints as well. This process can be continued until the entire panel surface is covered to the user’s satisfaction.

2.6. Adjusting the Projection Factor

The projection factor $p_f$ controls the distance between the two adjacent nodes (see Section 4), and is used to control the amount of gap or overlap between the viewpoints. Generally speaking, there should be no gaps between the viewpoints for complete coverage over the panel. However, for ease of understanding/preview, it is more helpful if there are some gaps between the outlines of the neighboring viewpoints. We set the projection factor to a default of $p_f=1.1$ for preview purposes and the user can adjust it by sending appropriate commands to change it in increments of 0.05.

Figure 10 shows the impact of adjusting projection factor. The default value of $p_f=1.1$ is shown in Figure 10a. On increasing it to $p_f=1.2$ in Figure 10b, we can see that the gaps between the neighboring viewpoint outlines have increased. Decreasing it to $p_f=1.0$ (Figure 10c) reduces the gaps to almost non-existent (in fact there will exactly be no gap for a developable surface) and decreasing it further to $p_f=0.9$ (Figure 10d) creates an overlap between the neighboring viewpoint outlines, which is favorable for complete coverage. Unfortunately, Figure 10d is difficult to interpret at first glance. The user is recommended to use the default value of $p_f=1.1$ and then decrease it based on the requirement to generate the preview.

After the coverage plan (viewpoints) have been generated to the user’s satisfaction, the travel path planning (Section 6) can be performed, and accordingly, the G-code or robot code for the manipulator to reach these viewpoints can be generated. Then the manipulator can be moved according to the code for doing the actual NDT inspection on the panel. These processes are not covered in detail in this paper.

3. Geodesic Paths on Surface Mesh

We now describe in detail the algorithm for finding the geodesic path of a given length L over the surface mesh given a starting point ${\mathbf{p}}_0$ and starting direction ${\widehat{t}}_0$. This section is already described in [10] but we repeat it here for more clarity. The length L corresponds to the movement of the camera from node to node as detailed in the next section and it is typically the distance to the next near-rectangle, in the length or width direction. This algorithm is taken from [11] with some modifications but it is worth repeating here and going into the details of how it is implemented. Briefly, it is an iterative algorithm that takes an initial point and direction, determines the next point and direction using a triangle intersection, and then repeats the process until the desired length is reached, or the boundary is reached, in which case it will extend the tangent over the edge up to the final length.

A pictorial explanation of the algorithm is shown in Figure 11. First, it is determined whether ${\mathbf{p}}_0$ lies inside a triangle, on an edge between two triangles, or on a vertex between several triangles. Accordingly, the tessellation normal $\hat{n}$ is computed. After that the direction vector ${\widehat{t}}_0$ is modified so as to make it perpendicular to $\hat{n}$ and then an intersection triangle is formed along the directions of ${\widehat{t}}_0$ and $\hat{n}$, with a length of d which is the geodesic path distance remaining to be covered (see Figure 11a). The size of the intersection triangle is also determined by ${\theta }_s$, which is the search angle. The search angle is used as a check to decide if the boundary is reached (rather than a continuation of the curvature) and is normally set to 30 degrees. Then we check for the intersection of this triangle with the triangles that ${\mathbf{p}}_0$ is a part of. There will only be one intersection line in all three cases (as we can see in Figure 11a–c), or no intersection at all in the case where there is no next triangle or the next triangle is folded at an angle greater than ${\theta }_s$ as we can see in Figure 11d. In both cases it means that the boundary has been reached, so we simply compute the final point with the distance remaining along the tangent vector at that point. The definition of ’boundary’ here is a sharp discontinuity in the triangle mesh as we can see in Figure 11d, and it is determined by the selection of ${\theta }_s$.

This algorithm also has two sub-functions: determining if p is a part of a triangle (including being on edges and vertices), and computing the tessellation normal n at that point. These two functions are described in Algorithms 1 and 2. Triangle–triangle intersection (computing the intersection line of two triangles) is described in [13] and not repeated here. The main algorithm is described in Algorithm 3.

Algorithm 1 is based on checking first whether the point is in the plane of the triangle, and then computing barycentric coordinates u,v to check whether it is inside the triangle or not. In the boundary cases where p lies on the edge or vertex, we would have at least one of $u=0,u=1,v=0,v=1$ (for vertices two of them would be true). Algorithm 2 computes the tessellation normal as per [11]. If the point is inside the triangle, the normal is the normal vector of that triangle. If the point is on an edge or vertex, the normal is the average of the normal vectors of all the triangles that this point p is a part of.

The inequalities mentioned in the above algorithms are usually evaluated in code with some tolerance built in, in order to account for numerical errors with floating point numbers (e.g., $u\le 1+ϵ$ rather than $u\le 1$). Also, for a closed surface not having a well-defined boundary such as a sphere or a torus, the geodesic path will simply wrap around the surface and continue. The starting point in Algorithm 3 can be any point on the surface. In practice, we start from the origin and project a ray upward along the z-axis. The intersection of this ray with the top of the mesh gives us our starting point.

Algorithm 1: Checking if point p is in triangle T

Algorithm 2: Calculating tessellation normal of point p

    The reason for using a geodesic path algorithm to generate a path of the surface is that geodesics are very well-defined curves that are easy to understand and model for parametric curved surfaces. However other techniques such as using slicing planes to generate a path could also be used and this is a topic that could be explored in the future. Also, the algorithm given in [11] is essentially an approximate way of computing a geodesic and this can also be improved. In the next section, we will see how the geodesic path algorithm will be used to generate the child nodes from parent nodes during node insertion.

An example output is shown in Figure 12b. A panel that is part of a launch vehicle propellant tank (it is an octant of a sphere with its top cut off) is shown in Figure 12a. We compute a set of geodesic paths on the 3D mesh model of this panel, starting from the origin point and length of 0.65 m, with directions of 5, 45, 75, 155, −75.

Algorithm 3: Generating Geodesic path sequence given input mesh, starting point, direction and desired path length

We note that the curve follows the surface, and it either completes within the mesh or upon reaching the boundary a last straight line is drawn off from it in the direction of the last direction vector up to the required length. In these cases, the end point of the geodesic path is outside the mesh.

4. Node Insertion Algorithm for Determining Viewpoints

We now detail our algorithm for inserting nodes, i.e., points on the surface from which viewpoints will be generated. These nodes are added in a tree-like fashion: first, the initial node is expanded in four directions then its child nodes are expanded, and so on until the node propagation stoppage criteria are met for all the nodes. This has been described in [10] but the stoppage criteria has been updated in this paper.

4.1. Nodes and Viewpoints

A new data structure is created for the node and its details are explained in

Table 1. Node data structure.

Data	Data Type	Explanation
p	vector(3)	Position vector of the node in R3 (on Mesh)
n	vector(3)	Outward normal vector of the node
t	vector(3)	Tangent vector of the node
b	vector(3)	Binormal vector of the node
i	integer	X-coordinate of the node in the 2D space
j	integer	Y-coordinate of the node in the 2D space
active	boolean	Whether the node is active or not

.

We note that the position of the node is recorded in both 3D and 2D space. In 3D space, the node’s position p is a point (x,y,z) on the mesh that is calculated by the geodesic path algorithm. In 2D space, the node’s position is not in units of distance but as integer (i,j) which denotes the projection from the initial point in the horizontal (tangent) and vertical (binormal) directions. The initial point is at (0,0) and expanding to the right creates a node at (1,0) and so on.

The normal, tangent, and binormal $\mathbf{n},\mathbf{t},\mathbf{b}$ are unit vectors that form a right-handed vector basis with its origin at the node (also called Frenet Serret reference frame), and the normal is oriented outwards from the mesh so that the viewpoint v can be located at some offset distance $d_0$ looking down on at the node. The view camera’s up-direction vector is oriented along the binormal direction. Figure 13 describes the node and view camera.

Figure 13 also shows us the projection of the viewing frustum on the surface. As we have assumed a quasi-flat approximation, we can say that the projection on the surface is approximately that of the viewing frustum. The horizontal projection $p_h$ and vertical projection $p_v$ can be thus calculated by the following formulas (Equation (1)).

$$\begin{array}{cc}\hfill p_h& =2p_f{d}_{0}tan({\theta }_{HFOV}/2)\hfill \\ \hfill p_v& =2p_f{d}_{0}tan({\theta }_{VFOV}/2)\hfill \end{array}$$

(1)

This can be understood by referring to Figure 14. We see in Figure 14b that on an ideal flat surface, the viewpoint can see a frustum of length $2d_{0}tan({\theta }_{HFOV}/2)$ in the horizontal direction (and $2d_{0}tan({\theta }_{VFOV}/2)$ in the vertical direction. Since we are using a quasi-flat approximation, we will assume that we need to move the same distance away to create another viewpoint with no overlap or gap. We add a projection factor $p_f$, which will decide whether the viewing areas have gaps in between them or overlaps. If $p_f>1$, there will be gaps (because we will be inserting nodes further away than their actual projection on the surface), and if $p_f<1$, there will be overlaps (Figure 14a shows gaps and Figure 14b shows overlaps). In practice, we want some overlaps between images but this makes the visualization of the coverage algorithm very confusing; hence, in this paper, we set $p_f=1.1$ in order to visualize each viewing area (corresponding to the node and viewpoint) separately on the surface.

The lengths $p_h$ and $p_v$ are important because these are our length inputs (L) to the geodesic path algorithm. We will find the child nodes from the parent node by moving by $p_h$ in a horizontal (tangent) direction and $p_v$ in a vertical (binormal) direction along the surface.

We also have two variables ${\theta }_{HFOV}$ and ${\theta }_{VFOV}$; these denote the horizontal and vertical field of view of the camera viewing frustum which we can see in Figure 13. These variables are not independent but related through the aspect ratio using Equation (2)

$${\theta }_{HFOV}={\theta }_{VFOV}\ast WIDTH/HEIGHT$$

(2)

where $WIDTH$ and $HEIGHT$ correspond to the sensor’s image size and are measured in pixels. For this paper, we use the parameters for an ELP IP camera with the camera’s parameters as per

Table 2. Camera  parameters.

Parameter	Value	Units
WIDTH	1280	Pixels
HEIGHT	738	Pixels
${\theta }_{VFOV}$	18	Degrees

. In this way, our virtual camera can be modeled for use in our algorithm, corresponding to a real camera. Note that focal length is not an independent parameter (it depends on sensor width, height, and field of view) and is not modeled here. The offset distance $d_0$ is a constant for all the viewpoints as we want all the images taken to be consistent and roughly the same size, and this distance is fixed at 30 cm for our algorithm. This gives us a viewing area of about 16.5 cm × 9.5 cm on the surface.

The nodes are inserted adjacent to one another until the node propagation stoppage criteria are reached for every node. The details of the algorithm are discussed in Section 4.3 and the stoppage criteria are described in the next section.

4.2. Node Propagation Stoppage Criteria

In this section, we will discuss the criteria for stopping the node propagation process. In our previous work [10], we observed that since the geodesic lines could draw close together or move apart on a highly curved surface like a hemisphere, too much overlap or separation of adjacent projections (bad behavior) could occur. The only stopping criteria at that point was to check if the newly placed node was outside the boundary of the mesh. In this paper, we introduce criteria for node stoppage to limit the view propagation so that the neighboring nodes do not converge or diverge from each other too much. We propose and explore two approaches here—the ’indirect’ approach which relies on normal vectors not diverging too much from the starting point, and the ’direct’ approach (which relies on managing the distance between neighboring viewpoints). Both the approaches are explained in the next section and the Results section shows the output from both methods so that we can compare them for many different use cases.

As we are not partitioning the surface, when we use a stopping criterion, only a portion of the surface may be covered before stopping. So we need to repeat the process from a yet-to-be-covered portion. Hence, in addition to the primary stopping criteria that we discussed in the previous paragraph, we need to also ensure that the current node propagation process does not move into an already existing set of nodes.

4.2.1. Indirect Method (Using Gauss Maps)

The idea behind this method comes from the concept of Gauss maps in differential geometry (Refer, e.g., [14]). The Gauss map is a map that maps any surface to a unit sphere, which transforms the unit normal at every point on the surface to the corresponding point on the unit sphere. If we draw the Gauss map for developable surfaces, i.e., plane, cylinder, and cone, we observe that the normals lie at the same point, along a great circle, or along a small circle, respectively (Figure 15a). On the other hand, for a highly non-developable surface such as a hemisphere, the Gauss map occupies a large area on the unit sphere. In fact, the Gauss map of a hemisphere is also a hemisphere (Figure 15b).

From the above observations, we define the concept of a near-developable patch on the surface as being an area where the normals diverge from each other by a maximum of some angle (say 60 degrees). In this region, we can still fit a planar grid onto the surface with no significant deviation. Thus, the stoppage criteria for the direct approach is arrived at. If the normal vector of the current node being inserted is more than a certain angle ${\theta }_N$ with respect to that of the starting node, we will stop the node propagation at this node.

The process is detailed in Algorithm 4.

Algorithm 4: Indirect method—Stoppage criteria

This method has a drawback that developable surfaces such as a cylinder will still stop the node propagation at that value of normal angle divergence ${\theta }_N$, even though the node propagation can still continue without any distortion (as the surface is developable). We will then need to take viewpoints from other starting points as well after saving the current set of viewpoints. This can be fixed to some extent by considering angular deviations in two directions, but as that also has some drawbacks, we do not explore it here. Instead, we use an alternative technique named the ’Direct method’, explained in the next section.

4.2.2. Direct Method

The idea behind this method is to check for convergence and divergence of the current node being inserted with respect to that of its neighboring nodes (if they already exist). In particular, if the node is located at (i,j) coordinates with respect to the start point (See Section 4, Table 1), we check the distance to its neighboring nodes at (i + 1,j), (i − 1,j), (i,j + 1), (i,j − 1) as seen in Figure 16. We compare the distance between current node (i,j) and nodes at (i + 1,j), (i − 1,j) with the horizontal projection on panel $p_H$ (Section 4, Figure 14) and the distance between current node (i,j) and nodes at (i,j + 1), (i,j − 1) to the vertical projection on panel $p_V$. If either of the two distances is less than or more than its prespecified limit, the propagation will be stopped.

The complete algorithm for implementing this criteria is given as follows (Algorithm 5).

Algorithm 5: Direct method - Stoppage criteria

4.2.3. Checking near Existing Saved Nodes

The final check for node propagation stoppage is to see whether the current node being inserted is close to any saved nodes or not. This is necessary to ensure that the node propagation is stopped if it ’runs into’ already existing nodes that were saved from a previous coverage planning process initiated from another start location.

This criteria is detailed in Algorithm 6.

Algorithm 6: Checking near saved nodes—Stoppage criteria

The value of $p_V$ is chosen conservatively as the minimum separation allowed between nodes, as there will always be an overlap if the node separation is less than $p_V$ (no matter which direction these viewpoints are angled, assuming $p_V<p_H$). Note that the Euclidean distance is used here with the assumption that the radius of curvature is large compared to the projection of viewing frustum on the surface.

4.3. Node Insertion Algorithm

We now detail our node insertion algorithm. A starting point is chosen on the mesh, its normal is determined and a direction is given as its starting tangent vector (This is given as an input by the user taking into consideration the panel shape). A set of nodes is created which is initially empty. We now insert the starting node into the set, with (i,j) as (0,0) and the binormal vector can be calculated from the normal and tangent vectors. Whenever we insert a node, it is initially set to active (An active node is one whose neighbors have not been fully expanded). We now continue adding nodes until there are no more active nodes. If an active node is found, we create new nodes along right, up, left, and down directions using the geodesic path algorithm, taking care to ensure that there isn’t already a node present at that position. The direction vector input depends on tangent or binormal vectors (for left/right or up/down) and the length input is $p_h$ or $p_v$ as described in Figure 14a.

We insert a new node at that point with the appropriate position, normal, tangent, binormal vector, and i,j set as active except in the case where the Node stoppage criteria are satisfied. in which case we make it inactive. After all the expansion possible is completed we set this node to inactive. We repeat this until all the nodes have become inactive which leads to a complete coverage of the area.

Node propagation stoppage criteria are explained in detail in Section 4.2 and is briefly recapped here. The criteria are that any one of the following are true:

• The newly inserted node lies outside the mesh;
• The newly inserted node is too close to an existing saved node;
• Either the ’Indirect method’ or ’Direct method’ is used to check the divergence of the nodes on the curved surface.

Node insertion in progress is shown in Figure 17. Currently, 12 nodes have been inserted so far. Figure 17a shows the position of the nodes on the surface of the mesh, the geodesic paths taken to get there, and also the corresponding viewpoints and intersection of viewing frustum with the mesh surface. Figure 17b shows the same node insertion progress in the 2D space of (i,j). The active nodes are shown in black and inactive nodes (which have been fully expanded) are shown in red.

The complete algorithm for node insertion is given in Algorithm 7.

Algorithm 7: Node Insertion Algorithm

5. Results

We run the Coverage Path Planning (CPP) node insertion algorithm for a variety of surfaces. Two of them are panels whose hardware is available to us from the Aerospace Industry (ISRO). One is of an octant shape (part of a propellant tank) and a cylindrical shape with cutouts, part of a launch vehicle structure. These two are listed as doubly curved convex and singly curved convex shapes, respectively. By rotating them upside down, we also obtain a doubly curved concave and singly curved concave shape, respectively. This is also a realistic use case in the industry as we need to do this kind of NDT inspection (thermography) for both sides of the panel. Three more panel shapes are considered here in mesh form (although we have no physical panels of that shape available to us)—A highly doubly curved convex surface (hemisphere), highly single curved in one direction and slightly curved in the other (nose cone shape), and a concave–convex structure (saddle shape). These examples are used to illustrate extensively the shapes of panels used in aerospace applications (see, e.g., Figure 18). Flat or near-flat shapes (e.g., airfoils, wings, empennage) and full cylinders (e.g., fuselage, engines) are not listed as they are basic variations of listed shapes. Note that all results here are produced with a default value of $p_f=1.1$. See Section 2.6 to see how the results change based on the projection factor.

The complete list of panel shapes is thus:

• Singly curved convex: Partial cylinder panel with cutouts, part of a launch vehicle.
• Singly curved concave: The same panel as singly curved convex but rotated 180 degrees about the Y-axis.
• Doubly curved convex: Octant shape, part of a propellant tank.
• Doubly curved concave: The same panel as doubly curved convex but rotated 180 degrees about the Y-axis.
• Doubly curved convex highly curved: Hemispherical shape.
• Nose cone: Highly curved in one direction and slightly curved in the other, with a singularity at the tip.
• Concave–convex: Saddle shape.

5.1. Singly Curved Convex Shape

This panel is part of a launch vehicle body. The cylindrical surface has 10,194 triangles and an upper surface area of 22,529 cm².

5.1.1. Indirect Method

We select two starting points using the indirect method to cover the entire surface.

5.1.2. Direct Method

The direct method can cover the surface with just one starting point.

We observe that both the indirect method (Figure 19) and direct method (Figure 20) give excellent results, although the indirect method requires two starting points owing to the limitation described in Section 4.2.1. The number of nodes generated in the indirect method is 172 compared to 161 using the direct method. This is because of overlap in the area covered by the first and second starting points as seen in Figure 19b.

5.2. Singly Curved Concave

5.2.1. Indirect Method

Two starting points were used to cover the surface, as seen in Figure 21.

5.2.2. Direct Method

As in the previous case, we observe that both the indirect method (Figure 21) and direct method (Figure 22) give excellent results, as this is the same panel but rotated to make it upside down. The indirect method requires two starting points just as in the previous case. The number of nodes generated in the indirect method here is 164—owing to a better selection of starting points, we have less overlap.

Given the viewing frustum area of 156.75 cm², we would expect 143 nodes for the cylindrical panel (singly curved concave/singly curved convex) surface coverage if the area was to be covered exactly. We see that we get more nodes than expected because many of these are near the boundary and cover very little of the surface. This also applies to the previous case of the singly curved convex surface as it is the same panel.

Both the singly curved convex and concave surfaces are equivalent to a curved plane; hence, there is no divergence of geodesic lines observed.

5.3. Doubly Curved Convex

This surface is an octant of a sphere (part of a propellant tank). It has 3046 triangles and an upper surface area of 13,644 cm².

5.3.1. Indirect Method

The indirect method does not cover the entire surface, leaving uncovered areas towards two of the corners. This is because of the limitation as described in Section 4.2.1. We omit multiple start points for full coverage here for simplicity.

5.3.2. Direct Method

We observe that both the indirect method (Figure 23) and direct method (Figure 24) give excellent results. Note that the geodesic lines start drawing close together (converging) towards the edges of the panel. This is due to the inherent doubly curved convex nature of the surface.

5.4. Doubly Curved Concave

5.4.1. Indirect Method

The indirect method requires two starting points for full coverage (Figure 25).

5.4.2. Direct Method

As in the previous case, we observe that both the indirect method (Figure 25) and direct method (Figure 26) give excellent results, as this is the same panel but rotated to make it upside down. The number of nodes generated in the indirect method here is 97, and in the direct method is 95, as compared to 87 nodes that would be expected for perfect coverage using the surface area of the octant-shaped panel. This is because of overlaps and wasted space at the boundaries.

We also note that in this case, the viewpoints diverge (or grow farther from each other) as we go from the center to the boundary. This is to be expected from the nature of geodesics on a doubly curved concave surface. It is in contrast to the converging nature observed in the case of the doubly curved convex surface.

5.5. Highly Double Curved—Hemispherical Surface

5.5.1. Indirect Method

For the coverage of the hemisphere using the indirect method, we need quite a few starting points. We have four main starting points (start points 1 to 4) corresponding to the four quadrants, then a fifth one for the top, then a sixth one to cover some gaps in between, and so on (See Figure 27). The complete coverage needed a total of 404 nodes using this method.

5.5.2. Direct Method

The advantage of the direct method in this case is that we can generate the entire coverage just in one click given a starting point (Figure 28a). In this case, what we observe is that the node propagation gets stopped at a certain region, but then it gets covered anyway by nodes coming from a different direction.

We also observe lots of overlaps (Figure 28) looking at the coverage from different angles by rotating the 3D view. These overlaps are caused when nodes from different directions of propagation bump into each other. Also, the notion of ’fitting a grid to the surface’ is essentially lost in this case because some distance after the start point, the nodes get propagated in ways that are not very intuitive. In the indirect method, by comparison, we can easily see how a grid is fitted to small patches of the sphere without much distortion. A total number of 421 nodes are generated using the direct method to cover the surface.

5.6. Nose Cone Surface

5.6.1. Indirect Method

Three starting points were used for coverage of the nose-cone shape (Figure 29).

5.6.2. Direct Method

We observe that both the ’Indirect method’ and ’Direct method’ perform poorly on the nose-cone surface (Figure 29 and Figure 30). It can be seen that geodesic lines ’wrap around’ the top point of the nose cone causing weird and non-intuitive node propagation. Moreover, the coverage area produced by the viewing frustum on the surface is no longer near-rectangular near the tip of the cone. In this case, effective coverage is hard to evaluate. The indirect method works better and more intuitively at least near the base of the cone. A total of 97 nodes were generated by the direct method as compared to 90 nodes by the indirect method for this surface.

5.7. Concave Convex Surface—Saddle Shape

5.7.1. Indirect Method

Two start points can be used to cover most of the saddle-shaped surface (Figure 31); however, some area near the top is left uncovered, which would need a few more start points to be fully covered (not shown here).

5.7.2. Direct Method

The direct method works oddly and non-intuitively for the saddle shape, as the convergence and divergence of geodesics along different directions causes the node propagation to leave large patches of the surface uncovered (Figure 32). A few more start points (not shown here) would be needed to cover these patches. In this case, also, we can say that the indirect method produces a more intuitive result.

5.8. Limitations of the Method

As we can see, the proposed method works reasonably well for all the panel shapes we are likely to encounter and the indirect method is more likely to give better results in the cases where the direct method is not working well. In general, for more complex shapes, the further we get from the central assumption that near rectangular patches are being fitted to cover the surface, the more likely that this technique will not work well; however, this would be true of any other algorithm as well. The other major limitation of this technique is that it is not fully automated but requires user input. The user has to select the starting criteria properly, which can take some experience.

5.9. Computational Analysis

It was observed that this algorithm is extremely fast and takes only a few seconds on a normal laptop. This means that the preview corresponding to a start point and direction chosen by the user can be generated nearly in real time without causing frustration to the user. We can see from Algorithm 7 that the computational complexity of the algorithm is $O\left(n^{2}M\right)$ where n is the number of nodes and M is the number of triangle faces in the mesh. This is because Algorithms 6 and 7 both iterate over the number of nodes (n) and Algorithm 3 iterates over the number of triangles (M), thus leading to the basic formula of $O\left(n^{2}M\right)$. Thus the algorithm runs in polynomial time. It can be noted that the number of nodes n can be computed approximately in advance by dividing the panel surface area by the surface area covered by each node. Hence, n depends on the camera parameters (field of view, distance from the surface, etc.) and also on the size of the panel.

While a direct comparison with optimization-based methods is outside the scope of this paper (due to the absence of literature addressing this specific topic of panel inspections), we can gain some insight into the speed efficiency of this method by looking at the performance of the drone-based coverage algorithm in [5] that is doing the coverage path planning over a church represented as a triangle mesh. In [5], the coverage planning algorithm takes 2.48 min to complete for the church model represented by just 200 triangles. In comparison, the method presented in this paper takes only a few seconds to cover the panels represented by over 10,000 triangles.

Further computational improvements can be made by using connectivity information between triangles (instead of searching the whole mesh every time in Algorithm 3) and also between nodes in the 2D parametric space. This would necessitate using a half-edge data structure for the mesh and additional changes. More computational improvements are possible by modifying Algorithms 6 and 7 to search only an appropriate subset of nodes instead of iterating over all nodes. However, because the current (‘inefficient’) algorithm still runs in just a few seconds on a normal laptop, the run speed was not considered a significant problem and the focus was more on simplicity rather than the efficiency of the algorithms. This can be explored in future work.

Since these results are obtained for models of actual panels used in the aerospace industry, we believe that our algorithm can directly be applied in this area for non-destructive testing of actual panels. The digital models of the camera and the surface are faithfully followed in order to get representative results relevant to the industry.

6. Travel Planning

The node insertion algorithm presented above generates the set of viewpoints, but it does not tell us how to move from one viewpoint to another in order to cover the surface. In general the problem of determining a minimum length motion path for a set of points such that each point is reached at least once is part of the Traveling Salesman Problem (TSP) which is NP hard. We should note here that in the type of application planned here (Thermographic NDT), the inspection time is much greater than the travel time, and thus there is no need to emphasize optimality in the path planning process.

The points generated from a single start point (using the indirect method as the stoppage criteria) are set out corresponding to a 2D grid. We propose a simple algorithm which would give a satisfactory solution to the problem. Note that this is only applicable to the case where these points are mapped roughly to a grid on the surface.

Figure 33 shows us the 2D space in which nodes have been mapped, each node representing a viewpoint. We proceed to develop a boustrophedon-like coverage path, starting from the right-lowermost node and moving up until all the nodes in that column have been covered. Then we move to the previous column’s uppermost node and change the direction of motion. This path planning is quite intuitive and so we do not describe the algorithm in more detail. Refs. [1,15] describe the boustrophedon cell decomposition methods which is what we follow but with a small difference. The ’cutout regions’ where there are no points in the 2D space do not correspond to ’obstacles’ in the conventional coverage path planning sense, since we can just pass over them (there is no physical obstacle). Hence we do not need to break down the region into cells and so on. We can see how the ’cutout region’ is handled in Figure 33, which shows the travel planning in process. In this figure, black dots denote the viewpoints (represented by their (i,j) values) already covered by the travel plan and red dots denote the viewpoints yet to be covered. Although this travel path is not necessarily optimal, it is good enough for practical purposes.

In the case where we have multiple start points, we can link the endpoint of the first travel path to the start point of the second travel path thus generated, and so on. Figure 34 illustrates this process.

For the case of path planning using the direct method for stoppage criteria, or to ensure optimality in general, we can implement the solution using the traveling salesman approach. Since solutions to the traveling salesman problem are well known in the literature, they are not explored here.

Since each node corresponds to a viewpoint, and each viewpoint can be represented by the position and orientation of the sensor, we can thus derive the motion of the machine or robot for moving to these viewpoints in sequence. We could thus generate CNC G-code or robot code for this purpose to be used in the actual application.

7. Conclusions

We have demonstrated a simple, easy-to-use, and fast algorithm for coverage path planning for non-destructive testing of panels used in aerospace applications using camera viewpoints. This uses a geometric method relying on generating nodes on the surface using geodesic path-finding, and it gives good results for representative sample panels used in the aerospace industry. The semi-automated approach relies on the user to set the start criteria, and two alternative approaches for setting node propagation stoppage criteria provide versatility and efficiency for covering the entire surface. This algorithm is much simpler than any optimization-based method and more suitable to be implemented in real-world applications.

Future work for this would include implementing the geodesic path-finding algorithm using a half-edge data structure for the mesh in order to make it more robust and computationally efficient. We could also calculate and demonstrate the actual coverage area on the panel for greater benefit to the user.

A further future scope of work would be to rigorously analyze the spatio-temporal complexity of the algorithm. An alternative method using an optimization and heuristic-based approach could also be implemented to solve the same problem of panel scanning. Metric comparatives could then be generated to qualitatively assess the performance of geometric methods and optimization methods.