Patent Application
20250014278
Aspects of the invention relate to:
Within the first of these a number of sub-sections cover the different elements up to and including BOLT (see below).
To convert a pointcloud (however derived) into an “intelligent” mesh where each separable planar surface (defined by break-lines) can be treated/processed individually ready for transfer (using a variety of industry standard formats) into other external software platforms eg Revit, Navisworks, Sketchup, etc. The process involves very considerable data compression thereby enabling the pointclouds to be operable on standard desk-top computers; where very large pointclouds are encountered they can be broken down using a methodology known as tiling (to be described).
Apart from the unique feature of producing an intelligent mesh (as explained separately) the Pointfuse processes retain all the relevant characteristics of the original pointcloud such that the accuracy of the final models relative to the original pointcloud can be represented on each section digitally or by using a heat map (by the application of standard deviation metrics).
The issues which are met in practice include
The solution ideally needs to have no gaps (so hole-filling is a fundamental requirement) and be “watertight” (to enable for example 3D printed models to be created from the Pointfuse output.
The algorithms created use advanced statistical sampling and probability theory specifically applied. Earlier attempts to produce a robust suite of software foundered because the known mathematical theory was applied to the individual points: the resulting code became unworkable to deal with the very large number of “special situations”. This application relies on applying the same underlying theory to whole planes rather than the points themselves.
As well as planes, the same underlying methodology can be applied to other surfaces and objects, for example, to cylinders and pipes.
The invention has numerous applications, including, for example, for use with Building Intelligence Modelling (BIM).
The invention provides solutions to the above and other problems of the prior art. Aspects of the invention are set out in the claims.
The features of the dependent claims may also be applied to other methods that create a representation of objects using a point cloud, for example, by deriving a surface mesh of planes or polygons, such as the method of WO 2014/132020 A 1.
In other words, the invention also provides a method of processing point cloud data of objects to create a representation of the objects, the method further comprising the features of any of the dependent claims, alone or in combination.
For example, the disclosure encompasses a method of processing a point cloud including point cloud data of objects to create a representation of the objects, the method comprising:
In this document:
A low curvature node that is adjacent to high curvature nodes along high curvature links and a t which the curvature does not change sign.
The line defined by the intersection of two distinct smooth surfaces.
A mesh link that belongs to only one polygon.
A node that belongs to an edge link.
A grid of cubes is superimposed over the region of space containing the point cloud. Any cube of this grid is called a grid cube.
The 3×3 cube comprising the primary cube and its 26 secondary cubes.
The vertexes of the polygon mesh.
A point in a point cloud.
Intersecting polygons. The links of the mesh correspond to one or more intersecting polygons.
Any individual grid cube.
A low curvature node that is adjacent to a boundary node along a contour.
One of the 26 grid cubes that are adjacent to the primary cube.
Non-overlapping cubes of a point cloud. A tile is 3D, though with airborne Lidar and some of the simplified figures here, you may see only 2D.
Either a vertex of grid cube or a vertex of the planar intersection polygon formed by the intersection of a surface and grid cube.
Pointfuse 1 converts a point cloud into planes (more strictly plane polygons) using the following steps, which are included in the recent Pointfuse patents (see WO 2014/132020 A1, the contents of which are incorporated herein by reference):
These steps result in a collection of plane polygons. At the borders of the point cloud, the polygon edges run along the grid cube faces. This looks fine if the plane is parallel to one of the cardinal planes, but otherwise the resulting polygons have unsightly zig-zag borders. In both cases, the polygon edges do not correspond with the position of the true point cloud border. Pointfuse 1 has an efficient statistically based method that corrects this error. This method was not included in the patent, and because it assumes the surfaces are planes, is not implemented in Pointfuse 2.
At this stage the surfaces are all still plane polygons. To display them on a computer screen, Pointfuse converts the polygons into a triangular mesh, using a standard triangulation method.
Finally, Pointfuse 1 can project the triangles of the 3D surfaces onto horizontal planes (plan view) or vertical planes (elevations). Pointfuse 1 implements Z-buffering to ensure that portions of (or all of) more distant triangles are correctly covered by nearer triangles.
Pointfuse 2020 consists of the following distinct but inter-operable components:
Pointfuse Kernel converts a point cloud into one or more separable meshes. It comprises the following steps.
The Pointfuse Tiler and Tilefuse enable Pointfuse to handle point clouds of unlimited size. The Tiler partitions the point cloud into one or more non-overlapping cubes called tiles. The Tiler then allocates those points into a collection of sub point clouds. Each sub point cloud contains the point cloud points that lie within the tile or lie within a narrow border region surrounding the tile.
The Pointfuse algorithm is applied individually (either sequentially or in parallel) to every sub point cloud, resulting in completely separate meshes in each tile. Each tile mesh may itself contain several separate surfaces. For example, one surface may represent the floor of a room, and another surface may represent a wall.
Tilefuse fuses the separate tile meshes into a single mesh in such a way that the separate surfaces within each tile are fused to the appropriate surfaces or surfaces in neighbouring tiles.
The Pointfuse Tiler partitions space into tiles. The length of each side of the tile is a multiple of the length of the grid cubes. By experimentation we have found that a tile size of 250 times the grid cube size results in fast run times. But other multiples, including non-integer multiples, can be used.
The Tiler allocates the point cloud points within each tile, and within a narrow overlap region surrounding the tile to a sub point cloud.
The overlap region extends a multiple of the grid cube length outwards from the tile in the positive and negative X, Y and Z directions. In one embodiment of the invention, the overlap region extended 5 grid cubes from the tile. Other multiples, including non-integer multiples, can be used.
Tilefuse is presented with a collection of tiles. Each tile contains a mesh which itself may be composed of sub meshes which are referred to here as surfaces. Each surface is distinct except that it may share nodes with adjacent surfaces in the same tile. Each surface may extend across the tile's six faces into neighbouring tiles. Tilefuse cuts each surface to ensure that the surface ends at the tile face. As a result of the surface cutting, some surface nodes now lie directly on the tile face.
Tilefuse now iterates through all the common tile faces (tile faces that belong to two tiles). Each common tile face has two adjacent tiles. Within each adjacent tile there is a mesh consisting of separate surfaces. If a surface has nodes or edges that lie on the tile face, it may be necessary to fuse those nodes and edges with the corresponding nodes and edges of a surface in the matching tile on the other side of the tile face.
Surfaces within each tile are cut by considering each triangle in turn. If all three triangle nodes lie inside the tile, or on the tile face, the triangle is left unchanged. If all three nodes lie strictly outside the tile, the triangle is removed from the surface. Otherwise, one of the triangle nodes must lie on the other side of the tile from the other two nodes. Two triangle edges must therefore cross the tile face. All such tile face crossing edges are identified. Usually they will belong to two triangles in the surface. All triangle edges that cross the tile face are split into two by adding a new split node at the position where the edge crosses the tile face. (Exception: if the position of a split node lies close to one of the edge nodes, that edge node is moved to the position of the split node and the edge is not split.) In
Once this process has been completed, the affected mesh polygons have either three, four or five nodes. If an affected polygon has only three nodes, no further action is required. Otherwise two of the polygon nodes are split nodes. If these two polygons do not already share an edge, the polygon is split along the line joining the two split nodes. The resulting two separate polygons lie on opposite sides of the tile face. The polygon that lies outside the tile is removed from the surface. If the remaining polygon has four nodes, it is split into two triangles, both of which lie inside the tile.
Once the external portions of the meshes have been removed, the meshes in adjacent tiles may be fused (that is combined) to form a single mesh. Meshes represent distinct surfaces so they are only fused if the angle between the planes of adjacent triangles are sufficiently small.
Nodes and edges to be fused are identified by establishing their proximity to sufficiently close nodes and edges in a matching surface on the other side of the tile face. For example
Let ak∈3 be the position of one or more nodes in a surface. And let b∈
3 and c∈
3 be the positions of the end nodes of an edge in the matching surface. The point on the straight line through b and c that is nearest to ak is:
If λk is sufficiently close to zero then λk is set to zero and ak will be fused with b.
If λk is sufficiently close to unity then λk is set to one and ak will fused with c.
Otherwise if 0<λk<1 then ak will be fused with a new node on the edge.
Otherwise ak will not be fused.
For example, in
The values λk and ak are stored in increasing order of λk for each edge b, c.
Each edge b, c that potentially can be fused lies along the tile face is the base of a triangle whose apex node d lies strictly in the tile or on the tile face. See
For example, in
Usually only single pairs of nodes are fused. Suppose the nodes being fused are A and B, with B being the node in the matching surface. Both nodes are moved to their common mean position. The matching node B is replaced by A in every triangle that B belongs to.
In the more general but rarer case, m nodes A1, . . . , Am in a surface and n matching nodes B1, . . . , Bn a matching surface are to be fused. In this case, all the nodes are moved to their common mean position and all nodes are replaced by A1 in every triangle that they belong to.
If the fused nodes belong to same surface (because the two surfaces were fused at an earlier stage) then no further action is required.
However if the fused nodes in the matching tile belong to a different then all the other nodes in the matching surface must also re-indexed to ensure that their labels are different from those in the current surface.
Consider a point cloud representing a surface that crosses a tile face. The portions of the surface on either side of the tile face are modelled by two completely independent processes in the two tiles. The two modelling processes separately generate two plane meshes that approximate the single point cloud surface. Because the two processes will in general work through the grid cubes in different orders, the two meshes will not be in perfect alignment.
Usually this misalignment is negligible. However when the surface is parallel, or nearly parallel, to the tile face, then even small discrepancies in the two meshes may result in the two surfaces crossing the tile face at significantly different locations.
To avoid this issue, meshes that are parallel or nearly parallel to the tile face, are first aligned to each other. This is done as follows:
Let the two meshes be labelled A and B.
Notice that, so far, analogous alignment changes have been made to both mesh A and mesh B. Steps 6 to 13 analyse mesh A, but make changes only to mesh B:
If a triangle PQR only has one edge PQ that is a cross link then the triangle is external if the node R lies outside tile B. In more detail, let the position of the triangle vertexes be p, q, r∈3, let u=r−p and let v be the external direction at P. Then triangle PQR is external if uTv>0. See
A grid of cubes is superimposed over the region of 3D space containing the point cloud. In what follows, a primary cube is any individual grid cube. Secondary cubes are the 26 grid cubes that are adjacent to the primary cube. The neighbourhood cube is the 3×3 cube comprising the primary cube and its 26 secondary cubes. See 3. Then the centroid of these points is
where the summation is over all point cloud points in the neighbourhood cube, and n is the number of point cloud points in the neighbourhood cube. The covariance matrix W∈3×3 of these points is:
where the superfix T indicates the transpose. The plane is computed in the form aT(x−3 is the plane normal and is the (unit) eigenvector corresponding to the smallest eigenvalue of the covariance matrix W. The plane can also be written in the form aTx=β where β=aT
Once a plane has been fitted to every grid cube, a vertex plane can be fitted at all the grid vertexes. (The vertex plane does not in general pass through the vertex, but is valid in the region surrounding it.) Each grid vertex belongs to eight grid cubes. In general a valid plane (called here a cube plane) has been fitted to each of these grid cubes. (There may be less than eight cube planes if some have been bound to be not valid). The vertex plane is an average of the valid grid planes and is computed as follows. Let the cubes planes be akTx=βk. Compute the symmetric positive semi-definite matrix A=ΣkakakT, where the summation is over the valid neighbouring cube planes. Let λ be the largest eigenvalue of A and let a∈3 be the corresponding (unit) eigenvector. Then the vertex plane is aTx=β where
where the summation is over valid the neighbouring cube planes. To see how this formula for β arises, note that by definition Aa=λa, that is
Similarly β=Σkwkβk.
A surface (treated locally as a plane) can intersect a cube at one or more of its vertexes and/or through one or more of its edges.
Consider the cube vertex Vi. Let the vertex plane at Vi be aiTx=βi. The displacement (or signed distance) of Vi from the vertex plane is di=aiTVi−βi. Then the surface is deemed to intersect the cube at Vi if |di|<0.001 h, where h is the resolution (the size of each grid cube).
Now consider the cube edge Eij joining cube vertexes Vi and Vj. If the surface intersects the cube at either Vi or Vj, it is deemed not to intersect the edge Eij. So in what follows it is assumed that the surface does not intersect the cube at either Vi or Vj. Then both di and dj are non-zero. The sign (the orientation) of the normals is arbitrary. For consistency it is therefore necessary to change the sign of the normal aj if aiTaj<0. This in turn changes the sign of βj and dj. We now treat the displacement as a linear function along the edge of Eij. If di and dj have the same sign (by construction neither can be zero) then the linear displacement function cannot be zero along the edge.
However if di and dj have opposite signs, then by linear interpolation the displacement is zero at the position V=(1−λ)Vi+λVj where λ=di/(di−dj).
The vertexes and edges that are intersected by the surface can be stored as a bit pattern as follows. Give the eight cube vertexes and 12 cube edges unique indexes in the range 0 to 19 as indicated in
If a cube grid does not possess a valid plane but has a vertex plane at each of its 8 vertexes, then those 8 vertex planes are used compute the intersection polygon, as described above.
The planes of cubes whose polygons do not possess a valid winding vertex set are recomputed as follows. Let the smoothing cube set consist of all cubes which do not possess a valid winding vertex set.
Any cubes which now have valid winding vertexes are removed from the smoothing cube set and the above steps 1 to 3 are repeated for each cube in the updated smoothing cube set.
We have not found it necessary to repeat steps 1 to 3 more than once.
This procedure is applied to the intersection polygons of pairs of adjacent cubes. It computes the distance between vertexes of the two polygons. The vertexes are deemed to be shared (the same) if the distance between them is less than a specified fuse tolerance. A suitable value for the tolerance is one tenth the length of a grid cube side. If the number of shared vertexes is less than three then no action is taken. Otherwise the two polygons must be coplanar. They are replaced by a single polygon that contains their combined vertexes (each shared vertex treated as one vertex) in the correct order around the combined polygon perimeter.
The edges of intersection polygons in adjacent cubes are compared. Polygons with more than one shared edge will have already been merged at the previous step. Any polygons without a shared edge are erased.
This trimming stage removes polygons that are both:
The smooth polygon creation stage ensures that surfaces meet exactly between adjacent faces of each grid cube. This step fuses together the nodes shared between the faces of grid cubes to make a contiguous surface.
Build change sets of polygon node indexes (sets of nodes that lie within the fusion tolerance of each other). All the nodes within a change set will be fused together and their indexes will be set to a single value. This includes nodes within the same polygon. It is possible that the change set will cause polygon butterflies (non-consecutive repeated node indexes).
Each node can fuse itself with any neighbouring nodes within the threshold distance. This step merges these neighbouring nodes together. This widens the fusion region but often resolves polygon butterfly issues.
Fuse the nodes in each change set together. Fusing nodes can cause the following issues:
This clean up step removes duplicate polygon nodes (short edges). This step can generate polygons with only one or two vertexes. These polygons are removed.
This clean up step removes duplicate polygons. That is polygons that consist of the same nodes, either in the same or reverse order.
The purpose of this calculation is to trim the polygon mesh to an estimate of the point cloud boundary.
This determines which polygons the active contour polyline can wander through to progressively move towards the point cloud boundary. An inward vector is calculated for each polygon in this work zone that determines which way the nodes of the active contour should travel. The work zone is built gradually, starting from the polygons connected to the outer contour of the mesh. Additional polygons are added to this set by gluing neighbouring polygons within 3 grid cubes distance of the outer contour.
The boundary edges of the mesh can be found by building a map of edges that link neighbouring polygons. The non-shared edges in this map identify the boundary edges. These unordered edges can be used to trace a path around the contour of the mesh to form a complete boundary polyline.
A map of nodes connecting unordered boundary edges must be built to trace the required edges around the contour. Outer contours of the mesh can legitimately touch at nodes where the mesh contains cavities. In cases where there are many possible outgoing edge links from a node, the edge whose connected polygon plane conforms most with the current surface normal is chosen.
Erode any polygons connected to the outer contour whose grid cubes do not contain any points. The outer contour polyline is subsequently updated after removing these polygons. This is performed to ensure the active contour can lock onto the nearby point cloud.
Polygon planes are computed from the vertices of the polygon. Each plane has a basis frame at the polygon centre with principal axes aligned with the inward vector, plane normal and a third unit vector, the binormal, orthogonal to both.
The work zone polygon planes and inward vector field are computed as follows:
This stage builds a local search horizon for each polygon within the work zone. Each polygon plane maintains links to its neighbours so that the nodes of the active contour polyline can freely move around the work zone. A 3-cube neighbourhood of polygon planes around the current one is gathered. Links are created to those polygons that are node or edge connected to the current one.
Each local horizon also stores the projected point cloud points for the active contour to lock onto. The horizon of projected points is built as follows.
The boundary of the point cloud must be determined in a different way for these polygons to prevent trimming polyline interaction and unintentional erosion. The active contour polyline is disabled for edges connected to bridge polygons. A polygon in the work zone is labelled as a bridge type if: all of its neighbouring polygons either touch the outer contour (forms a single polygon strip), or it is part of a two-adjacent work zone configuration whose polygons are connected to the outer contour and have opposing inward vectors (2-polygon side by side strip).
The active contour polyline is a 3D polyline that rides over the work zone surface in the inward direction towards the boundary of the point cloud. The algorithm ensures that the polyline is always contained within the work zone with its nodes projected onto the surface. The polyline successfully traverses both flat and curved regions of the mesh by iteratively smoothing the active polyline and driving it towards the boundary of the point cloud.
It is an iterative technique that progressively increases the fraction that the active contour moves towards the point cloud boundary. Only 5 iterations are required to refine the contour to the boundary shape.
For each polyline node get its previous and next nodes. In the general case the smoothed node position is moved to the mid-point of the line formed between the previous-current edge midpoint and the current-next edge midpoint. The endpoints remain fixed for open polylines that do not form a complete loop. The number of polyline nodes remains the same after the smoothing operation.
The nodes of the active contour polyline must lie within the work zone at all times. In the general case the polyline node moves between polygon planes of the work zone. The horizon of polygon planes associated with the node at its last valid position is used to search for its current position.
The smoothing step can also pull the active contour polyline outside of the mesh itself especially in concave regions. In such circumstances a step is performed to backtrack the movement of a polyline node until it lies within the work zone again. The horizon of polygon planes is searched until the line intersects an edge of its plane polygon.
Each active contour edge moves independently.
The inward distances to the point cloud boundary are estimated from both of the edge's endpoints. Different inward distances at the edge endpoints allow the active edge to turn towards the point cloud boundary.
The active contour nodes separate when determining the boundary.
The node positions are averaged later to recover a connected polyline.
The polyline edge's inward vector is determined from the cross product of the polygon plane normal and edge direction. The edge's inward vector is transformed into the frame of the polygon plane (on 2D plane) and conditioned to lie on the same side as the polygon plane's inward vector.
A frame is created at each endpoint of the active edge to probe the horizon's point cloud samples (inward vector, projected node position and complement direction).
The horizon of point cloud samples is transformed into the frame of the active edge endpoint. The near and far distances of samples along the edge's inward vector are computed.
The point cloud boundary is computed by finding the distance from the far baseline that accounts for 90% of included points. A binary search strategy is used to bracket the threshold distance that gives the appropriate number of included points. See
The drive target point for the active edge's endpoint is then easily determined from the boundary distance.
In one embodiment a caching mechanism is used to prevent too much re-computation of the point cloud boundary whenever the active edge does not move too much.
The target can be pulled outside of the mesh itself, especially in concave regions. In such circumstances a step is performed to backtrack the movement between target and current node position until it lies within the work zone again. The horizon of polygon planes is searched until the line intersects an edge of its plane polygon.
At the start of the procedure the active contour polyline has edges that are aligned with the faces of the grid cubes. The first stages of the algorithm are biased towards polyline smoothing to remove the serrations seen in the outer contour. This quickly aligns the polyline edges with the point cloud boundary. It makes sense to increase the fraction the contour drives towards the boundary target as the number of iterations increase.
Short active contour edges tend to crumple and turn too sharply whenever the polyline is forcibly constrained to the work zone boundary. A minimum edge length constraint is employed to prevent serious issues like this. At the end of each iteration any short edges are removed, with those edges connected to it having their endpoints moved to the removed edge's mid-point.
Extend open polylines to the work zone edges to ensure that the outside area of the trimmed surface is removed correctly.
Any polygon in the work zone should only be trimmed by a single polyline. The part of this trimming polyline that resides in the neighbourhood can be generated in three steps.
Iterate over the edges of polygons in the work zone. Gather any cutting polylines from either polygon connected to the edge and determine the points where the trimming lines pass over the edge. Each polygon edge has an associated basis frame that is used as a common space to perform intersection tests.
The candidate cutting points are stored in the edge data. Close candidate edge intersections are then merged.
Nodes are then added into the polygons that connect to the edge.
This stage orders the newly added intersections with the polygon edges along the trimming polyline. This synchronization makes it easy to replace parts of the polygon's edge sequence.
Traverse the ‘trimming edge’ polyline based on its endpoint Tags to build the edges of the polygon that need replacing:
Convert the line segments back to a polyline format that has no duplicate nodes. The vertical positions of any points that are added into the interior of the polygon are interpolated in the plane normal direction between new polygon edge points.
This stage splits the polygon into two parts by finding where the trimming polyline crosses the contour. The inside part of the polygon containing the point cloud is combined with the part of the trimming polyline that crosses the polygon between entry and exit points.
The intersection of two distinct smooth surfaces defines a line which is called here a break line. Note that the smooth surfaces are not necessarily planes and therefore break lines are not in general straight lines.
The previous work flow stages have:
Each polygon resides within a single grid cube, called the polygon's primary cube. The mesh is called smoothed because each polygon has been fitted through the point cloud points contained in the 27 grid cubes centred on the polygon's primary cube. As a result the mesh gives an excellent fit to smooth surfaces such as planes, pipes or cones. However a much less satisfactory fit is obtained at break lines, for example on stairs, or at the corners of buildings. The purpose of the break line calculation is to improve the fit at such locations by identifying the intersection between the multiple surfaces and modelling the surfaces separately.
Let pi∈3 be the position of a mesh node. A mesh node pa∈
3 is said to be adjacent to pi if the two nodes are consecutive vertexes of a mesh polygon. Let
i be the set of all mesh nodes that are adjacent to pi. The covariance matrix at pi is defined as:
A surface trihedron consisting of three orthonormal vectors ui, vi, wi∈3 is computed at each node pi. The vector wi is computed as the unit eigenvector that corresponds to the smallest eigenvalue of wi. The other two vectors ui and vi are the remaining two eigenvectors of Wi.
If we treat the mesh nodes as lying in a smooth (that is, continuously differentiable) surface, then wi can be interpreted as the surface normal at pi and the other two vectors ui and vi form a basis of the tangent plane to the surface at pi. See
The set i of nodes pa that that are adjacent to pi is also called the depth 1 set of nodes connected to pi and is written
i(1). More generally,
i(r), the depth r set of nodes connected to pi, is the set of nodes that belong to
i(r-1) or are adjacent to any node in
i(r-1).
Each mesh polygon's vertexes lie on the edges of its primary cube, and the polygon's edges lie across the cube's faces. Therefore all the polygon edges are parallel to one or other of the cardinal planes, and therefore are orthogonal to one or two of the coordinate axes. All polygon edges that are orthogonal to the X axis are called “X edges”, all polygon edges orthogonal to the Y axis are “Y edges”, and all polygon edges that are orthogonal to the Z axis are “Z edges”. Note a polygon edge that lies along a grid cube edge will be simultaneously orthogonal to two of the coordinate axes.
In
“X polylines” are polylines consisting entirely of X edges. “Y polylines” and “Z polylines” are defined analogously. They are surface contour lines corresponding to constant X, Y and Z values.
Note that a node will always lie on at least two contour lines, and sometimes on three such lines.
If the pi∈3 is adjacent along X edges to two other nodes pa∈
3 and pb∈
3, then the normal curvature at pi (in the plane through pi that is orthogonal to the X-axis) is computed by the formula κX=rTwi/rTr where wi∈
3 is the unit surface normal at pi and r=p0−pi, where p0∈
3 is the position of the centre of the circle fitted through the three nodes pa, pb, pi.
If pi is adjacent along an X edge to only one node pa, then the normal curvature is computed by the formula κX=2 (pa−pi)Twi/(pa−pi)T(pa−pi).
The normal curvatures Ky and Kz are computed analogously.
A node will always have at least two normal curvatures.
The curvature tensor at pi is estimated as follows. Let ui, vi, wi be the surface trihedron at pi. Every node pk∈(3) in the depth 3 set of nodes connected to pi, can be written as
where xk=ΔpkTuk, yk=pkTvk, zk=ΔpkTwk, and where Δpk=pk−p0. Compute the unit vector:
This corresponds to a unit vector in the tangent plane. Compute the curvature associated with this direction as
More generally, if the curvature tensor has coefficients
in this coordinate system, then the normal curvature along the unit direction ({circumflex over (x)}, ŷ) is given by:
It is natural to compute these coefficients by linear least squares. That is we minimize
The first step of the break line calculation is to locate nodes with high surface curvature. Break lines are likely to occur close to such nodes.
Because the calculation of the curvature tensor is relatively expensive, the normal curvatures are first estimated using the appropriate formulae of Section 4.3.5. If the absolute values of any of the normal curvatures is greater than the tolerance κTOL, then the curvature tensor is computed as explained in Section 4.3.6. The node pi is considered to have high curvature if either of the absolute values of the eigenvalues of the curvature tensor
are greater than κTOL.
The curvature value of every high curvature node is stored in a high curvature node table. The table contains a separate curvature value for each contour that passes through the node and the links associated with that curvature. Such links are called high curvature links.
A suitable choice for the numerical value of curvature tolerance κTOL is essential. If κTOL is set too high then subsequent steps in the break line calculation may create break lines where none exist. Similarly, giving κTOL too small a value will cause some break lines to remain unrecognized.
A good compromise is obtained by setting κTOL=2/9 h where h is the size of the grid cubes (the resolution).
Isolated nodes are defined as nodes that do not themselves have high curvature, but are connected on both sides along a contour by high curvature nodes. The isolated nodes are added to the high curvature node set, together with the links that connect them to high curvature nodes. See
Consider a contour line passing in order through three adjacent nodes p0, p1 and p2 on a surface. Let
and let û and {circumflex over (v)} be the corresponding unit vectors. Compute the vector product
In detail, let û=(ux, uy, uz), {circumflex over (v)}=(ux, vy, uz) and w=(wx, wy, wz). Then
But û and {circumflex over (v)} lie along the same contour line and therefore either are both X edges, both Y edges or both Z edges. If they are both X edges, then their x-components uk and uy are both zero, in which case both wy and wz are zero, and the sign of the curvature at p1 is defined as the sign of wx. Similarly if û and {circumflex over (v)} are both Y edges then wx and wz are both zero and the sign of the curvature at p1 is defined as the sign of wy. And if û and {circumflex over (v)} are both Z edges then the sign of the curvature at p1 is defined as the sign of wz. Notice that the sign of curvature, calculated in this way, does not depend on the orientation of the surface normal.
The sign of curvature is calculated for every high curvature node (including isolated nodes). A separate sign of curvature is computed for each contour passing through the node.
If the node is a point of inflection along a contour (that is if along the contour the curvature has different signs on either side of the node) then the curvature value corresponding to that contour is removed from the high node curvature table entry for that node.
Boundary nodes are low curvature nodes that:
Boundary nodes are added to the high curvature table, together with their high curvature links.
A double boundary node is one that is adjacent to another boundary node along a high curvature link. Double boundary nodes are removed from the high curvature table, together with the link connecting them.
The surface normals at all double boundary points are recalculated so that they are orthogonal to the line joining them. In detail, let pa, pb∈3 be a pair of boundary nodes and let their corresponding surface normals be wa and wb. Then these surface normals are recomputed to become
where u=pb−pa and αa=−uTwa/uTu and αb=−uTwb/uTu.
An edge link is a mesh link that belongs to only one polygon. An edge node is a node that belongs to an edge link.
The previous node table is used to trace a chain of nodes along a contour. It associates each node in the chain with its predecessor nodes along the contour. Because there are three possible contours through each node, the previous node table may associate nodes with three different previous nodes.
A root node is a low curvature node that is adjacent to a boundary node along a contour. The previous node table is seeded by allocating storage space (but initially without associated previous nodes) for all root nodes, double boundary nodes and curvature change nodes (nodes at which the curvature changes sign, that is the nodes on opposite sides along the contour have different signs). Exception: any nodes that are also edge nodes are not added to the previous node table.
The previous node table is grown by the following process. Identify all high curvature nodes that are both not in the table and adjacent, along high curvature links, to nodes that are already in the table. Once all of the adjacent nodes have been identified, each adjacent node is added to the previous node table in two places. In detail, suppose node A is already in the table and node B is an adjacent node that has high curvature along the link joining A and B. Then:
This process is repeated until no further suitable adjacent nodes remain.
The previous node table implicitly defines chains of consecutive nodes along a contour. Apart from the two end nodes of the chain, each chain node has two predecessor nodes. The two end nodes, having only one processor node, are readily identified as root nodes. The chain root table lists the root node pairs associated with each chain node. There is one root pair for each chain that the chain node belongs to. The chain root table is compiled as follows. Starting at a chain node, which by definition has two predecessor nodes, the chain is traced in both directions until the two root nodes are reached at the end.
Each chain lies along a contour and therefore in a plane orthogonal to one of the coordinate axes. The surface normal at a root node is called skew if it is nearly orthogonal to the coordinate axis. More precisely, the root node is skew if êrTŵa>0.9, where êr∈3 is the direction of the coordinate axis orthogonal to the chain and ŵa∈
3 is the surface normal at the root node.
If both root nodes are skew, the root nodes and all the chain nodes between them are removed from the chain root node table.
If only one root node is skew, an attempt to find a better replacement root node pair and chain is made as follows. Suppose that the current chain nodes listed in order along the chain are a, c0, c1, c2, . . . , cn, b where a and b are the root nodes and a is skew. If one of the chain nodes c, also belongs to another chain, then the current chain is shortened to become cr, . . . , b and the root nodes of the chain become cr and b. The chain root node table is updated accordingly. If more than one of the current chain nodes also belongs to another chain, then the one nearest to b (along the chain) is used to replace a.
The position and surface normal of each root node define a plane through that root node. Each chain has two root nodes and therefore has two associated root planes. The key idea of the break line calculation is to move each chain node (each node between the two root nodes) by projecting it onto the nearest of the two root planes. The surface normals at each moved chain node is also changed to equal the normal of the corresponding root plane. See
Now consider two chains A and B. Suppose that node R is a root node of chain A and is a chain node (not a root node) of chain B. Then node R is called a dependent root node. Its position and surface normal (and therefore its plane) will be projected onto the nearest of two root nodes of chain B. It is essential that the projection of node R occurs before R is used as a root node, so that the resulting chain node positions are consistent. To ensure this, chain nodes are only projected if both their root nodes are not dependent. If any of the newly projected chain nodes are root nodes of another chain, these root nodes (which by definition are dependent) are flagged as having been resolved. Once both root nodes of a chain are either not dependent or have been resolved, the chain nodes can be projected onto the nearest root plane. See
Once the chain nodes have been projected, every chain node lies in one of the two root planes. A break link is any link in the chain whose ends do not both lie in the same root plane. If there is only one break link chain (which is usually the case) a break node is inserted in the break link, splitting it into two links. The two end nodes of the link being broken are called side nodes. The break node is positioned at the intersection of the two root planes. (Exception: if one of the nodes of the break link is sufficiently close to both root planes, it becomes a break node and is moved to the intersection of the root planes. In this case the break link is not split. In this case one of the side nodes is identical to the break node.)
In
An affected node is a node that is not in a chain but is adjacent to at least two chain nodes that have been moved. Each of these adjacent nodes will have been projected to a root plane. The affected node is projected onto the nearest root plane.
In
Clone split nodes are break nodes that are sufficiently close to each other (for example within 0.2 times the grid cube size)) and have a common side node. They usually occur at the intersections of three distinct surfaces. (See
A binary is a node that only has two adjacent nodes. If a binary node is
Every break node is associated with the two root planes. The edge direction is the unit normal that is orthogonal to both plane normals. See
Polygons containing exactly two break nodes are split into two separate polygons along the line joining the two break nodes provided the break nodes are not adjacent and the angle between their edge directions is less than 30 degrees. Note that the two break nodes are contained in both resulting polygons and are adjacent in those polygons. In
A polygon that contains exactly one break node is called a break polygon. Two break polygons form a break polygon pair if they share at least two nodes but not break nodes. That is, between them, they contain two break nodes.
A break polygon pair is isolated if its two constituent polygons are not members of other break polygon pairs.
For example in
An isolated break polygon pair is split by replacing the two constituent polygons by a single combined polygon, and then splitting the combined polygon along the line joining the two break nodes.
Isolated break polygon pairs are split in two passes:
For example in
In both passes, the polygons are only split if the angle between the two break edges is less than 30 degrees and if the angle between the line joining the two break nodes and each break edge is less than 45 degrees.
A break polygon triplet are three polygons A, B and C where:
For example in
Break polygon triplets are split by joining the three polygons together to form a single composite polygon, and then splitting the composite polygon along the line through its two break nodes.
In this context a graph is a topological data structure consisting of nodes and links between pairs of nodes. A break line node is a graph in which all the nodes are break nodes.
An end node is a node that belongs to only one link. Since the links represent break lines, end nodes should only occur at locations where the point cloud ends. Otherwise they represent a location where the break line should be extended. One situation in which this can occur is that of triple nodes, where three break lines need to be extended to meet at a corner.
An end polygon is a polygon that contains one end node of the break line graph, but no other break nodes. A triple node T is a node that is not a break node but belongs to (is a node of) three distinct end polygons. Each of the three end nodes is associated with two root planes, making a theoretical total of six planes in all. However some of these root planes may be identical. If the system of six linear equations corresponding to the root planes has rank 3, then a unique intersection position p0∈3 can be computed by linear least squares. If the position of the triple node T is sufficiently close to the intersection position p0, then:
For example in
A polygon that contains more than two break nodes is called a corner polygon. Corner polygon clusters are sets of corner polygons that share at least one node (not necessarily a break node) with another corner polygon. However a corner polygon cluster will often contain only one polygon.
In
Each individual corner polygon P of a cluster may be adjacent to (share at least two nodes with) a polygon A that:
The set of all such adjacent break node polygons is identified for the cluster. Once the entire set has been identified, it is added to the cluster to form an augmented cluster. Thus only a single layer of adjacent polygons is added, and the cluster does not grow indefinitely.
In
Next the polygons adjacent to the augmented cluster are investigated. Each edge of such an adjacent polygon is called common if the edge is shared with a polygon of the augmented cluster, and is called separate otherwise. The adjacent polygon is said to be embedded in the cluster if the sum of the lengths of its common edges is greater than the sum of the lengths of its separate edges. The set of all embedded polygons is identified. Once the embedded set has been identified its polygons are added to the augmented cluster to form a set called a nearly convex corner cluster.
By construction, the nearly convex corner cluster will contain at least three break nodes. Each break node will be associated with two root planes, making a theoretical total of six root planes. If the rank of these root planes is three, they are solved by linear least squares to give a unique intersection point. This step is identical to that used in splitting triple nodes.
If the computed intersection point is sufficiently close to at least one of the nodes of the nearly convex corner cluster, then the intersection point is treated as a new break node and all polygons of the convex cluster are replaced by a set of new corner polygons, where each new corner polygon contains the new break node plus two existing break nodes. The two straight line segments joining the new break node to the two existing break nodes are both edges of the new corner polygon. The remaining edges trace a connected path from one existing break node to the other break node.
When this step is applied, the mesh consists of polygons, some of which may be non-planar. Polygons with more than three nodes are split into triangles in such a way as to minimize the sum of the cosines between each triangle normal.
At the start of the surface simplification step, the mesh consists of a very large number of triangles, many of which are smaller than a grid cube. Although the mesh is an excellent representation of the point cloud surface being modelled, the mesh requires so much memory to store it, that at least for larger models, computers struggle to handle it. This is particularly so if the mesh is passed to third-party software which may also be using the available memory for other purposes. The surface simplification step therefore attempts to reduce the number of triangles.
In what follows, an edge is called a break edge if:
An edge is that is not a break edge is called a smooth edge.
A triangle is called skinny if any of its angles is smaller than a given tolerance. In one embodiment of our method a tolerance of 5 degrees was used.
The surface simplification process consists of the following steps:
Each distinct triangle mesh is simplified by collapsing certain edges. That is, the two nodes that comprise the edge are replaced by a single node positioned near to the mid-point of the edge. The new node replaces the existing nodes in all the connected edges (the mesh edges that contain one of the nodes being replaced).
Every mesh edge is considered as a possible candidate for edge collapsing. Three properties are computed for each edge for use in the mesh simplification. The three properties are:
The collapse position is the location of the single node which would replace the edge if it were collapsed. The ‘maximum height from plane’ is a measure of how far the collapse position is from the planes of the surrounding triangles. (The collapse position and the maximum height from plane are defined in detail below. In particular, in some circumstances the maximum height from plane computation may be deemed invalid.) Because the positions of an edge's end points may change as the mesh simplification progresses, an edge's maximum height from plane and collapse position may be recomputed several times. Each time they are successfully recomputed, the edge's job number is incremented by one.
Edges whose ‘maximum height from plane’ computation is valid and whose maximum height from plane value is less than a given threshold are added to a map of edges to be collapsed. The threshold is called the planar patch fitting tolerance. A suitable value is one fifth of the grid cube size. The map key is the ‘maximum height from the plane’. The other entries are the candidate edge (defined by its end nodes) and the edge's current job number.
Each edge in the map of edges is considered in turn, starting with the edge having smallest maximum height from plane.
If the edge no longer exists in the mesh, or if the job number stored in the map does not equal the edge's job number, the entry is deleted from the map. The latter happens if the edge's maximum height from plane has been recomputed more recently than the current map entry.
Otherwise the edge's maximum height from plane and collapse position are recomputed (because the position of the edge's two nodes may have changed). If the computation is invalid the entry is deleted from the map. If the edge's new maximum height from plane is valid and remains below the planar fitting tolerance, the edge is collapsed and replaced by a node at its collapse position, and the current entry is deleted from the map. The ‘maximum height from plane’ values are recomputed for all the edges connected to the new node. If the recomputed value is valid and is below the planar fitting tolerance, the connected edge is inserted into the map of edges to be collapsed. Note that this means that an edge can occur in more than one place in the map, however the job number identifies the most recent one.
Once an entry in the map has been processed it is deleted from the map. The edge collapse process is continued until the map of edges to be collapsed is empty.
This section explains how the maximum height from plane is computed for an edge AB. As part of the process, the position of the collapse position is also calculated. A mesh triangle is called a connected triangle if its nodes contain at least one of the end nodes A and B. Let the plane of each connected triangle Δk be written as akX=βk, where ak is the plane unit normal, βk is the plane scalar constant, and the suffix k is an integer identifying the triangle. Let be the set of identifiers of all triangles connected to AB. Let
be the identifiers of all triangles connected to edge AB that contain only one of the end nodes A and B. Note that
is a subset of
. For each triangle Δr in
, compute a shape position pr∈
3 as follows. Suppose that triangle Δr has nodes P, Q and R, where R is one of A and B. Form the equilateral triangle PQS where the new node S lies on the same side of PQ as R. The shape position is given by the location of S. See
An edge of a triangle is called connected if it shares one or two nodes with AB. If one of a connected triangle's edges is AB, then all three of its edges are connected. Otherwise the triangle has two connected edges.
A push down list of planes is constructed as follows. The individual entries of the push down list are stored as pairs (ak, βk). Each connected triangle is visited in turn. The connected triangle's plane is pushed once for each of the triangle's connected edges that is also a break edge or a boundary edge (an edge that belongs to only one triangle of the surface). Note that the push down list will often contain the same plane more than one. This is because:
A second push-down list of planes is also constructed. The planes in the list are called end stop planes. Each connected triangle is considered in turn. Within each connected triangle, each edge is considered in turn. If the edge is connected to AB (that is, at least one of the edge's nodes are A or B) and if the edge is a break edge or a boundary edge (an edge that belongs to only one triangle) then the triangle plane is pushed onto the list of end stop planes. The purpose of the end stop planes list is to prevent a break edge or boundary edge being moved (at least not far) if the edge is collapsed.
The collapse position x∈3 is computed so as to minimize the sum of squares
where in the second summation on the right is over all entries in the push down list of end stop planes, k(r) is the identifier of the rth entry in the list, m∈3 is the mid-point of the nodes A and B, and w∈
is a weighting factor. In one embodiment of the invention the weighting factor is set to 0.01. The collapse position balances the conflicting requirements that it lies close to:
In one embodiment, the minimum of the objective function V is given by the solution of the linear equation Ax=b, where the matrix A∈3×3 and the vector b∈
3 are computed as
where n is the number of entries in . If A is not of full rank. the maximum height from plane computation is flagged as invalid.
To prevent the edge collapse from moving connected triangles on top of each other, the maximum height from plane computation is flagged as invalid if the collapse position does not lie on the same side of the base of each connected triangle as that triangle's shape position.
The maximum height from plane is the maximum distance of the collapse position from the planes of all the connected triangles.
Break edge runs are consecutive break (or boundary) edges. They stop at open edges or forks. An open (ended) edge is one that is not connected to any other break (or boundary) edge. The seed triangles are triangles that contain an open edge. (By construction each seed triangle can only contain one open edge). Each open edge therefore has two seed triangles.
Each seed triangle is grown into a nearly planar surface. The surface consists of a contiguous set of triangles. The triangles are contiguous in the sense that each triangle in the set shares an edge with another triangle in the set (unless the set consists of only one triangle.) The planar surface is grown using a push down list of edges. Edges are pushed onto the list and then popped (removed) after they have been inspected. The nodes of every triangle within the nearly planar surface lie within a tolerance of the seed plane (the plane of the seed triangle).
The growing process is as follows:
Step 2 is repeated until the frontier set is empty.
By construction, the nearly planar surface only contains one break or boundary edge. This is the seed edge, the open edge that forms one edge of the seed triangle. The nearly planar surface has its own border edges, that is edges that belong to only one of the triangles in the nearly planar surface. One of these is the seed edge. The border edges of the nearly planar surface are any edges that belong to only one triangle of the surface. (The term border is used to distinguish it from boundary edges of the mesh.) Border edges are traced from each open end of the seed edge, to construct a contiguous chain of border edges, until the chain reaches an end node of another break edge run. Because this chain of border edges may have branches, it is necessary to trace all possible edges, forming a tree of border edges. Note that each as end node has two trees, one for each seed triangle. Both the root node and the leaf nodes of the tree are end nodes. Each root node and leaf node are connected by a contiguous sequence of border edges, none of which are break or boundary edges. The distance between every root node and leaf node pair is the sum of the length of the chain of border edges. The end node pair with shortest distance is connected by adding every border node in their chain to the set of break edges.
Surface texturing is concerned with colouring the surfaces of the mesh surface to inform the point cloud being modelled. Note that here colouring includes RGB and intensity data that may be stored for each point in the point cloud. It may also include derived point information such as the distance of point from the mesh.
The colour information is constructed as a 2D surface called a texture. As well as its (x,y,z) 3D position coordinates, each vertex in the mesh has 2D (u, v) texture coordinates. There is one pair of texture coordinates for each mesh triangle that the vertex belongs to. The texture coordinates define the vertex's position within a bitmap. The colour, intensity or other texture value of any point on the surface of the mesh triangle is computed by linear interpolation in the bitmap.
The bitmap itself is constructed in two main steps. First the triangles of each surface are unwrapped onto a 2D surface. Secondly, the colour or other value of each pixel within the bitmap is accumulated by projecting the point data from the point cloud points associated with each triangle. Both processes are described below.
The process for optimally arranging the triangles within a rectangular bitmap is well known and is not described here.
Note that Pointfuse meshes have been constructed as one or more distinct surfaces. The unwrapping process partitions the triangle meshes of each surface into distinct subsets called cells. At the end of the unwrapping process, each triangle will have been allocated to one of these cells. The wrapped versions of each triangle in a given cell all lie in a single plane and do not overlap each other. In general, if the surface has high curvature, it will be unwrapped into more than one cell.
Each mesh node has two distinct related positions: its wrapped (x, y, z) position in 3D space and its unwrapped (u, v) position in the cell plane. The cell centroid is the mean of the unwrapped positions of the nodes in the cell.
In what follows, triangles are called neighbours if they share a common edge. The unwrapping process is:
where α=uT(c−q) and β=vT(c−q). That is, p is the projection of q onto the cell plane.
If the wrapped position p lies outside all triangles in the current cell, then:
Otherwise remove the neighbour triangle from and go to Step 4.
The unwrapped surface is coloured by projecting the point cloud points onto the interior of the original (wrapped) mesh triangles. The corresponding pixel is then coloured in the unwrapped version of the triangle. Two complications arise:
These are dealt with as follows:
In simple terms, the colour accumulator sums the colour values of every point within a non-empty pixel and then divides that sum by the number of points in that pixel to give the average pixel colour value. In practice, the bitmap image is mixed by distributing the colour information in a small region (the “influence window”) surrounding each non-empty pixel, and computing the average of these distributed values.
In one embodiment of this method, the influence window consists of the 9 pixels surrounding the central non-empty pixel. A separate triangular weighting function is constructed in the X and Y directions. Each function has its maximum value (one) at the projected position of the point cloud point, and has a width of two pixel lengths. The weight used in each pixel is the value of the weight function at the centre of that pixel. A combined weight is now calculated in each pixel in the influence window by multiplying X and Y weights together. Finally, the combined weights are normalized by dividing them by their sum.
The non-empty pixels are used to colour the empty pixels. In one embodiment, the process makes use of a filled set which initially consists of all non-empty (and therefore coloured) pixels. A frontier set, consisting of all the empty pixels that are neighbours of coloured pixels, is constructed. All the pixels in the frontier set are given the average of the colours of their neighbouring filled pixels. Once all the frontier pixels have been coloured, they are moved to the filled set and new frontier set is constructed. The process is repeated until any remaining empty pixels do not have filled neighbours.
Pointfuse Space Creator converts Pointfuse meshes into floor plans of user-selected storeys of a building. For the purposes of this note, each floor plan is treated as a simple graph consisting of nodes and edges. Each edge connects two nodes and represents a wall. Edges are called adjacent if they share a node and are called disjoint otherwise. Each floor plan can be exported for use by third party software. If the angle between two adjacent edges is approximately but not exactly 90 or 180 degrees, the two edges may sometimes be interpreted as disjoint by the third party software. To prevent this from happening, the positions of the wall plan nodes are tweaked, before export, so that all approximate right angles become exact right angles and all approximate straight angles become exact straight angles.
A connected node triplet Tijk consists of three nodes {i, j, k} such that node i is connected to both node j and node k. Let pi∈2 be the position of node i, let uij=pj−pi and let
The two unit vectors ûij and ûik define the directions of two walls that meet at node i.
The node angle θijk (in radians) between the two walls is defined by
The node angle is considered to be approximately a right angle if
and approximately a straight angle if
where θTOL is the target angle tolerance in radians. We have found that a suitable target angle tolerance is 5 degrees, that is θTOL=π/36 radians.
The problem of tweaking the positions of the wall plan nodes so that node angles that lie within θTOL of their target angles can be treated as the constrained optimization problem
where pi(0) is the starting position of node i and N is the number of nodes and where the constraints are
where τijk is the target angle (either π/2 or π) of the connected node triplet Tijk. Only connected node triplets whose angles lie within θTOL of one of the target angles are included in the constraints.
The solution algorithm is described in terms of the constrained optimization problem
subject to m nonlinear constraints
c(x)=0
where x∈n and x(0)∈
n are n-vectors and where c:
n→
m The problem variables x and the position vectors pi are related as follows. Let pi=(Xi, Yi) and pi(0)=(Xi(0), Yi(0). Similarly write x=(x1 . . . xn) and x(0)=(x1(0) . . . xn(0))T. Then n=2N and for i∈{1 . . . N}
Let x(0)∈n be a given n-vector. This section describes a fast (locally quadratically convergent) numerical method for finding a point in a nonlinear implicit manifold that is locally closest to x(0). The implicit manifold
⊂
n is defined by the m nonlinear equations
where c: n→
m are a system of m nonlinear continuously differentiable functions. That is
A locally closest point in is defined as a local solution x*∈
n of the nonlinearly constrained optimization problem
where the objective function
is the sum of squares of the distances of each variable xi from its starting point xi(0).
The special structure of this optimization problem allows it to be solved without explicit reference to second derivatives or Lagrange multipliers.
The solution method is iterative. Let x(k) (k≥0) be the starting point of the kth iteration. (In one embodiment of the method, the zeroth iteration starts at the given point x(0), but this is not essential.) The idea of the method is to linearize the constraints (1) about x(k) and to set the next trial solution x(k+1) equal to the unique point on the linearized constraints that is nearest to x(0).
In detail, the constraints linearized about x (k) are
where c(k)=c(x(k))∈m is the m-vector of constraint functions values evaluated at x(k) and Ak=∇c(x(k))T∈
m×n is the Jacobian matrix of the constraint functions evaluated at x(k) and Δx(k)=x−x(k)∈
n. It is assumed that at least one of the entries of Ak is non-zero.
Because the linearized constraints (3) may not be consistent we pre-multiply (3) by AkT to obtain the normal equations
By the QR factorization theorem there exists an orthogonal matrix Qk∈m×m and a permutation matrix Pk∈
n×n such that
where Uk∈r×r is upper triangular and invertible and upper triangular, Mk∈
r×(n−r) and r is the rank of Ak. If Ak is zero then r=0 and the method cannot be used. If Ak is of full rank (r=m) then (5) is
Equations (6), (7) and (8) below remain valid in this special case except that the zero sub-matrices and v(k) are omitted. The equations after (8) are unchanged. Substituting (5) into the normal equations (4) gives
Pre-multiplying (6) by Pk−T=Pk gives
where Δy(k)∈r and Δz(k)∈
n−r are related by
and b(k)∈r and v(k)∈
m−r are related by
Multiplying out (7) gives
The top row of (9) can be simplified as
Note that the bottom row of (9)
provides no further information because it is (10) pre-multiplied by MkT. Equation (10) can be written as
where y−y(k)=Δy(k) and z−z(k)=Δz(k) Equation (11) explicitly partitions the problem variables x into dependent variables y and independent variables z. Equation (11) can be further re-written as
Pre-multiplying (12) by Uk−1 gives
Note that because Uk is upper triangular the term s(k)≡Uk−1 b(k) in (13) can be computed by solving Uks(k)=b(k) by back substitution rather than by explicitly inverting Uk. Similarly the matrix Wk can be computed by solving UkWk=Mk. In this case, each of the n−r columns of Wk are computed separately; the jth column wk(j)∈r of Wk is computed by solving Ukwk(j)=mk(j) where mk(j)∈
r is the jth column of Mk.
Rearranging (13) gives
Now the objective function (2) can be written
Using (14) to eliminate the dependent variables y from (15) gives
Multiplying out the first term on the left hand side of (16) and rearranging gives
is the unit matrix of order n−r. Note that Gk is strictly positive definite and is therefore invertible. A person skilled in the art will recognize that Gk is the reduced Hessian matrix of the linearly constrained optimization problem being solved: minimize the sum of squares objective function ƒ(x) subject to the normal equations (4).
The gradient vector of ƒ(z) with respect to the (n−r)-vector of independent variables z is
Setting the gradient vector to zero, dividing by 2 and rearranging gives
Of course someone skilled in the art will compute z by solving the linear equations (18) numerically (for example by Gaussian elimination or Cholesky factorization) rather than by explicitly inverting Gk.
To obtain the corresponding values of the dependent variables y, one substitutes z into the linearized constraints (10). Subtracting z(k) from equation (19) gives
Substituting Δz(k) into (10) and rearranging gives
Because Uk is upper triangular, the value of Δy(k) can be obtained from (20) by backward substitution.
The new trial solution x(k+1) is given by
Let (x, R)={u∈
n: ∥u−x∥2<R} be the open ball with centre at x∈
n and R>0 is the radius. Here ∥ ∥2 is the Euclidean or L2 norm. An n-vector x*∈
n is a local solution of (P1) if c(x*)=0 and there exists a positive R>such that ƒ(x*)≤ƒ(x) for all x∈
(x*, R).
Our experience is that the sequence of trial solutions x(k) always converges rapidly to a local solution x*∈n of (P1), probably because the starting point x(0) is sufficiently close to x*, however there is no theoretical guarantee that convergence will occur unless additional measures are taken.
For example, in one embodiment, the simple update equation (21) is replaced by
where 0<ak≤1 is called the step length and is chosen in such a way as to ensure convergence.
In detail, let g(x) be the sum of the squares of the constraint function values at x. That is
where g(k)=g(x(k)). But by construction
showing that unless x(k) already satisfies the constraints, Δx(k) is a descent direction of g at x(k). In particular, provided only that Ak is non-zero, it follows that if x(k+1) is computed using (22) then
for all sufficiently small step lengths αk>0.
Initially set αk to one, and progressively multiply αk by a constant reduction factor 0<ω<1 (for example ω=0.5 or ω=0.1) until (23) is satisfied.
It remains to explain how the derivatives of the constraint functions are calculated. The position vector of node i is
The constraint associated with connected node triplet Ty can be written as
where τijk is the corresponding target angle and
It follows that
As in WO 2014/132020 A1, aspects of the invention can be carried out using suitable devices and apparatus, such as a scanning module, processing module, point cloud database, or computer etc.
Elements of the description can be combined with each other, and with elements of WO 2014/132020 A1.
| Number | Date | Country | Kind |
|---|---|---|---|
| 2108778.8 | Jun 2021 | GB | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/GB2022/051555 | 6/17/2022 | WO |
