CA3108547A1 - Collaborative 3d mapping and surface registration - Google Patents

Abstract

Subject matter regards generating a 3D point cloud and registering the 3D
point cloud to the surface of the Earth (sometimes called -geo-locating"). A method can include capturing, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of the object, registering the overlapping subsections to each other, and geo-locating the registered overlapping subsections.

Description

CLAIM OF PRIORITY
[0001] This patent application claims the benefit of U.S. Provisional Application Serial No. 62/975,016, filed February 11, 2020, which is incorporated by reference herein in its entirety.
TECHNICAL FIELD

[0002] Some embodiments described herein generally relate to generating a three-dimensional (3D) mapping and registering the generated 3D mapping to a surface.
BACKGROUND

[0003] Generating a 3D point cloud can be resource intensive. The generation of the point cloud can include gathering two-dimensional (2D) images (e.g., satellite imagery, ground imagery (images taken from a camera on the ground), or an elevation therebetween) and performing photogrammetry, performing a light detection and ranging (LIDAR) scan, a human generating a computer-aided design (CAD) drawing, or the like. Each of these techniques has associated difficulties, tolerances, and resource requirements, which can be cost or resource prohibitive.
BRIEF DESCRIPTION OF THE DRAWINGS

[0004] In the drawings, which are not necessarily drawn to scale, like numerals may describe similar components in different views. Like numerals having different letter suffixes may represent different instances of similar components. Some embodiments are illustrated by way of example, and not limitation, in the figures of the accompanying drawings.

[0005] FIG. 1 illustrates, by way of example, a conceptual block diagram of an embodiment of a technique for 3D point cloud registration, such as with error propagation.

[0006] FIG. 2 illustrates, by way of example, a diagram of an embodiment of a technique of performing the operation of the processing circuitry.

[0007] FIG. 3 illustrates, by way of example, a conceptual block diagram of an embodiment of a collaborative system for 3D point cloud generation.

Date Recue/Date Received 2021-02-10

[0008] FIG. 4 illustrates, by way of example, a diagram of an embodiment of a system for generating the registered 3D point cloud.

[0009] FIG. 5 illustrates, by way of example, a diagram of an embodiment of a volume that is being mapped by UAVs.

[0010] FIG. 6 illustrates, by way of example, a diagram of an embodiment of a system for 3D point cloud generation and geo-registration.

[0011] FIG. 7 illustrates, by way of example, a diagram of embodiments of operations that can aid in geo-locating a 3D point cloud.

[0012] FIG. 8 illustrates, by way of an example, an embodiment of a system 800 for 3D point set registration and merging.

[0013] FIG. 9 illustrates an example diagram of an embodiment of the relationship between ground point coordinate estimates V3 and corresponding 3D
data set observations V,.

[0014] FIG. 10 illustrates an example of an embodiment of a bundle adjustment operation.

[0015] FIG. 11 illustrates, by way of example, a diagram of an embodiment of a method for 3D point set generation and registration.

[0016] FIG. 12 illustrates, by way of example, a block diagram of an embodiment of a machine in the example form of a computer system within which instructions, for causing the machine to perform any one or more of the methodologies discussed herein, may be executed.
DETAILED DESCRIPTION

[0017] Aspects of embodiments regard improving point cloud generation or registration of the point cloud to a geo-location of the Earth or other surface.

[0018] FIG. 1 illustrates, by way of example, a conceptual block diagram of an embodiment of a technique for 3D point cloud registration, such as with error propagation. A system of FIG. 1 includes input 102 provided to processing circuitry 114. The processing circuitry 114 generates output in the form of a 3D point cloud 128.
The input 102 can include one or more of an image or video 104, a light detection and ranging (LIDAR) scan 106, a three-dimensional (3D) computer-aided drafting (CAD) model 108, satellite imagery 110, or other data 112. The input 102 can be processed into two or more 3D point clouds before being provided to the processing circuitry 114.

Date Recue/Date Received 2021-02-10

[0019] The image, video 104 can include a red, green, blue (RGB), infrared (IR), black and white, grayscale, or other intensity image. The image, video 104 can include a video that comprises frames. Photogrammetry can be performed on the data of image, video 104, such as to generate one of the 3D point clouds.

[0020] Photogrammetry can include performing a geometric bundle adjustment on the two-dimensional (2D) images, to register the geometry of the 2D images to each other. The bundle adjustment can adjust geometry of an image of the images to be consistent with geometry of other images of the images. The geometry can be defined in metadata, such as by using rational polynomial coefficients (RPC). Other image registration techniques are possible. For example, 2D images not previously associated with the 3D point cloud can be registered to the 3D point cloud. Tie points can be identified between each 2D image and the 3D point cloud. The geometry of each image can be adjusted to match the 3D point cloud by using an affine transformation.

[0021] The LIDAR scan 104 can be generated by illuminating a target with a laser light and measuring the reflected light with a sensor. Differences in laser return times and wavelengths can then be used to make one or more of the 3D point clouds.
This is because these differences can be used to determine distance to the object off of which the light was reflected and returned to the sensor. LIDAR 3D point clouds often have no intensity information. For the LIDAR case, and others, it can be useful to attribute the 3D point cloud with intensity or color data from an image that covers the same area of interest. Further discussion on this point is provided elsewhere herein.

[0022] The CAD model 108 is a human-designed (e.g., with or without computer aid) 2D or 3D model, such as the blueprint of a building. The CAD
model 108 is defined by geometrical parameters and readily adjustable by a human using a computer.

[0023] The satellite imagery 110 includes images generated at high altitudes from one or more cameras of a satellite or satellites. The satellite image provides a nadir, or near-nadir view of a geographical location. The nadir view is of a point on a celestial sphere directly below an observer (e.g., the point on the sphere closest the observer).

[0024] Other sources 112 can include other man-made measurements (e.g., thermal imaging), computer aided measurements, or other data that can be used to generate a 3D point cloud of the geographical region.

Date Recue/Date Received 2021-02-10

[0025] The processing circuitry 114 can include hardware, software, firmware, or a combination thereof, configured to implement the operations of registering 3D
point clouds to each other. Hardware can include one or more electric or electronic components configured to electrical signals to indicate results of the operation. Electric or electronic components can include one or more transistors, resistors, capacitors, diodes, inductors, switches, logic gates (e.g., AND, OR, XOR, negate, buffer, or the like), multiplexers, power supplies, regulators, analog to digital or digital to analog converters, amplifiers, processors (e.g., application specific integrated circuits (ASIC), field programmable gate array (FPGA), graphics processing units (GPUs), central processing units (CPUs), or the like), or the like electrically coupled to perform the operations.

[0026] The operations for 3D point cloud registration can include 2D/3D
point cloud conversion and normalization at operation 116, 3D point cloud geo -registration at operation 118, and 3D point cloud fusion and adaptive filtering at operation 120. These operations are discussed in more detail below.

[0027] The operation 118, in general, can include determining a scale 122 factor adjustment, a rotation 124 adjustment, and a translation 126 adjustment between the 3D
point clouds to be registered. The operation 118 can include using an iterative, normalized cross-covariance technique, that minimizes a least squares difference between tie-points ground control points (GCPs), or the like. This is discussed in more detail below. The result of the registration is a registered 3D point cloud 128 that inherits the best errors (smallest errors) in the 3D point cloud inputs.
Again, more detail is provided below.

[0028] FIG. 2 illustrates, by way of example, a diagram of an embodiment of a technique of performing the operation 118 of the processing circuitry 114. As previously mentioned, more detail can be found below. In general, tie points 212, 214, 216, 218 between the various 3D point cloud inputs can be identified. In the example of FIG. 2, the 3D point clouds include the LIDAR scan 106, the CAD model 108, the image, video 104, and the satellite imagery 110. More or fewer 3D point clouds can be used. The tie points 212, 214, 216, 218 are data points that correspond to a same geographic location. For example, a corner of a structure, a high elevation point, or the like, can make for a good tie point. The tie points 212, 214, 216, 218 can be used to determine how to adjust the corresponding 3D point clouds to be registered to generate Date Recue/Date Received 2021-02-10 the 3D point cloud 128. The tie point 220 corresponds to the registered location of the tie points 212, 214, 216, 218.

[0029] In performing the operation 118 an accounting of the registered points 222, and the corresponding data that was used to the generate the registered points 222 can be performed. The registered points 222 can include an x value 224, a y value 226, a z value 228, metadata 230, and source pointers 232. The metadata 230 can include source reference vectors 234 indicating data source (e.g., the image, video 104, LIDAR
scan 106, CAD model 108, or satellite imagery 110, etc.) from which the registered points 222 were determined and the source data 236 of those sources.

[0030] A system, device, or method to implement the techniques of FIGS. 1 or 2 can provide a hybrid, multi-source, multi-resolution 3D point cloud creation and enrichment through multi-source, multi-modal 3D point cloud ingest and fusion.
The techniques of FIGS. 1 or 2 can provide an ability to create a more comprehensive (as compared to prior techniques), multi-resolution, 3D hybrid point cloud by combining a sparse 3D point cloud with a dense 3D point cloud, a low-resolution 3D point cloud with a high-resolution 3D point cloud or the like. The techniques of FIGS. 1 or 2 can provide an ability to fill in a missing section of a 3D point cloud, such as by using another 3D point cloud.

[0031] The techniques of FIGS. 1 and 2 can provide an ability to replace or augment noisy/low quality 3D point cloud sections with high fidelity data from other location-relevant 2D and 3D data sources. The techniques of FIGS. 1 and 2 can provide an ability to detect error and correct a 3D point cloud through multi-3D point cloud cross-correlation & validation. A resulting hybrid 3D point cloud can preserve, through metadata, data source lineage that allows users to leverage metadata (e.g., pixel color and intensity, object classifications and dimensions) from fused sources.

[0032] The techniques of FIGS. 1 and 2 provide a hybrid 3D point cloud-derived location intelligence (e.g., detection, localization, and classification of permanent or stationary objects in a scene or environmental conditions (e.g., power or phone lines, radiation areas that may not be detectable by the onboard sensors of an unmanned aerial vehicle (UAV) or an unmanned ground vehicle (UGV)) that can aid in obstacle avoidance and planning of complex mapping operations. The techniques of FIGS.
1 and 2 can provide an ability to incorporate fusion criteria, such as trustworthiness of the data source or classification levels.
Date Recue/Date Received 2021-02-10

[0033] The techniques of FIGS. 1 and 2 provide a user-controllable filtering and criteria, such as classification level of a data source or trustworthiness of a data source, quality, or age. The techniques of FIGS. 1 and 2 provide an ability to control which areas (e.g., rooms within a building) and/or objects get included or excluded in a resulting, hybrid 3D point cloud. The techniques of the FIGS. provide an ability to link other measurements (e.g., temperature, radiation, noise level, humidity) to each data point in the resulting hybrid 3D point cloud (e.g., via a graph, multi-dimensional array, hash table, dictionary, or linked list representation).

[0034] The extensible, hybrid cloud source vector 222 allows the resulting hybrid 3D point cloud to also store non-sensitive (e.g., unclassified) and sensitive (e.g., classified) point cloud data in the same cloud. Hybrid point cloud access controls, for example, can be enforced by controlling the visibility of the source data vectors 222, or encrypting the sensitive source data elements 236 of the hybrid point cloud.

[0035] FIG. 3 illustrates, by way of example, a conceptual block diagram of an embodiment of a system 300 for 3D point cloud generation. The 3D point cloud generated using the system 300 can be used as an input 102 to the processing circuitry 114, for example. The system 300 as illustrated includes unmanned aerial vehicles (UAVs) 330, 332, 334, 336, with imaging devices (indicated by diverging lines 338).
The imaging devices can include intensity or non-intensity imaging devices.
Examples of intensity imaging device include an RGB, grayscale, black and white, infrared, or another camera. Examples of non-intensity imaging devices include LIDAR, or the like.
The UAVs 330, 332, 334, 336 can be programmed to capture image data of an object 340 or a geographical region of interest. The UAVs 330, 332, 334, 336 can cooperatively capture sufficient image data of the object 340 to generate a 3D
point cloud of the object 340.

[0036] The concepts presented herein regarding UAVs are not limited to drone-/UAV-collected data. Concepts can also apply to UGVs, unmanned vessels, and manual-/human-driven data and collection methods.

[0037] FIG. 4 illustrates, by way of example, a diagram of an embodiment of a system 400 for generating the registered 3D point cloud 128. The system 400 as illustrated includes 3D mapping scheduler 440, a collaborative 3D object mapping operation 442, 3D point clouds 444 from the UAVs 330, 332, 334, 336, a 3D
point cloud database (DB) 446, the processing circuitry 114, and the registered 3D
point cloud 128 (both from FIG. 1).

Date Recue/Date Received 2021-02-10

[0038] The 3D mapping scheduler 440 can command the UAVs 330, 332, 334, 336 what tasks to perform and when. The 3D mapping scheduler 440 can change the task or time to perform the task after a mission has begun. The UAVs 330, 332, 334, 336 can communicate and alter the task or timing on their own. The UAVs 330, 332, 334, 336 can be autonomous or semi-autonomous.

[0039] The 3D mapping scheduler 440 can provide the UAVs 330, 332, 334, 336 with a task that includes a geographic region to be modelled, a resolution of the model to be generated, a technique to be used in generating the model (e.g., color image, satellite image, radiation or temperature scan, LIDAR, etc.), or one or more time constraints in performing the tasks. The UAVs 330, 332, 334, 336 can operate to satisfy the schedule and constraints provided by the scheduler 440. While illustrated as a centralized unit, the scheduler 440 does not need to be a centralized unit.
The scheduler can be distributed across multiple resources (e.g., UAVs) and run locally/on-board (e.g., as an agent) to provide a distributed dynamic scheduler. Such an implementation can include the UAVs 330, 332, 334, 336 communicating and planning tasks among themselves.

[0040] These operations of the UAVs 330, 332, 334, 336 are part of collaborative 3D object mapping operation 442. The result of the 3D mapping operation 442 can be a 3D point cloud 444 from the UAVs. As previously discussed, the discussion regards UAVs, but also applies to other manned or unmanned vehicles, such as ground, water, air vehicles, or a combination thereof.

[0041] The 3D point cloud 444 can be stored in a point cloud database 446. The point cloud database 446 can include a memory device for storing one or more point clouds, such as the 3D point clouds 444 or the registered point cloud 128. The 3D point clouds 444 from the UAVs 330, 332, 334, 336 can be provided to the processing circuitry 114. The processing circuitry 114 can register the point clouds 444 to generate the registered point cloud 128.

[0042] FIG. 5 illustrates, by way of example, a diagram of an embodiment of a volume that is being mapped by UAVs, such as three of the UAVs 330, 332, 334, 336.
Currently, part of the volume, as indicated at 502, is partially mapped. The different patterns on voxels or subsections of the volume indicate which UAV has mapped the subsection (if a UAV has mapped the subsection). The scheduler 440 or the UAVs 330, 332, 334, 336 can determine which UAV will map the currently unimaged subsections based on a variety of criteria (e.g., object priorities, speed, cost, distances, available Date Recue/Date Received 2021-02-10 sensors, required scan resolution, remaining flight times/battery life, UAV
health status, etc.). An example assignment of the mapping is provided at 504.

[0043] An advantage of one or more embodiments of FIGS. 1-5 can include one or more of: enabling collaborative, concurrent, inside-out and outside-in 3D
mapping of objects of interests (e.g., buildings, tunnels, or other objects) using a swarm of two or more autonomous, semi-autonomous, or man-controlled drones; providing an ability to leverage pre-existing 3D point clouds and 3D point clouds generated through other means (e.g., CAD models, photogrammetry) to speed up mapping process, preventing replication of effort, increasing fault-tolerance, and/or reducing cost;
generating a fused 3D point cloud that can be composed of multi-resolution 3D point clouds;
providing a collaborative mapping and hybrid 3D point cloud creation process supports dynamic filtering based on user-definable criteria (e.g., exclude basement from mapping, specified resolution, specified time frame, specified technique, or the like);
efficiently delegating mapping objectives / areas of interests using pre-assigned mission planner or assigned at run time by a dynamic scheduler; providing a dynamic scheduler that enables ad-hoc re-tasking of mapping objectives in response to changing mapping objectives, environmental conditions, or drone(s) health status (e.g., battery status, sensor or motor failure, data storage constraints, change in weather conditions, or the like); providing a swarming solution with autonomous mapping drones, since drones can negotiate / self-organize and potentially trade their 3D mapping assignments/objectives based on a variety of criteria or conditions (e.g., tasking from a failing drone gets automatically re-assigned to another qualified drone);
change in mapping priorities/order; auto-selection of closest drone(s) to targeted area(s); auto-selection of drone with fastest or most qualified sensor payload; selection based on drone's health status (e.g., remaining battery life or storage, strongest downlink signal, or the like); providing an ability to incorporate 3D no-fly zones and exclusion area definitions in mission planning and scheduling; providing an ability to explicitly task drones or allow autonomous, onboard decision making (e.g., based on mission objectives and/or current conditions and location) to adaptively switch between sensors (e.g., Lidar, IR, or high definition (HD) cameras), turn sensors on and off, and adjust sampling rates in flight to speed up mapping, optimize storage utilization, speed up onboard analytics, extend flight time/battery life, or reduce electronic signature;
providing an ability to link in UTM (Unmanned Traffic Management) systems (e.g., Federal Aviation Administration (FAA), Airmap) in 3D mapping planner, dynamic Date Recue/Date Received 2021-02-10 scheduler, and platform's real-time flight control system to avoid air space violations and collisions with other aircraft and unmanned systems; providing an ability to 3D map in large area or complex object in multiple stages (e.g., first use a survey drone to rapidly generate a coarse, low-res 3D point cloud of the area of interest, then use a geo-spatial portioning technique to subdivide the resulting 3D model (e.g., by floors, volume, openings (windows, doors), voxels, or sectors), and then task the survey drone(s) to explore certain areas / sub sections (sequentially or in parallel) within the coarse/quick scan-generated 3D model; or providing a 3D divide-and-conquer approach, enabling parallel, multi-drone and priority-driven 3D exploration and mapping.

[0044] FIG. 6 illustrates, by way of example, a diagram of an embodiment of a system 600 for 3D point cloud generation and geo-registration. The system 600 as illustrated includes an operator 602, the UAV 330, and a geographical region 604 to be mapped. The operator 602 can operate the UAV 330, or the UAV 330 can operate autonomously or semi-autonomously to generate image data of the geographical region 604. In some instances, the UAV 330 can have a global positioning system (GPS) or the like that informs the UAV 330 of its location relative to the surface of the Earth. In these instances, the GPS coordinates can be used to register the 3D data to the surface of the Earth. In other instances, however, the UAV 330 does not have such a system.
Note it is not necessary to have GPS data for the entire flight. GPS data for a portion of the flight can be used to register data from the entire flight. When the UAV does not have the positioning system, some other techniques can be used to register the 3D point cloud generated using the image data from the UAV 330 to the surface of the Earth.
Operations of the other techniques are provided regarding a method 700 provided in FIG. 7. Note that not all operations of FIG. 7 are required or even useful in all situations.

[0045] FIG. 7 illustrates, by way of example, a diagram of an embodiment of a technique 700 for registration of a 3D point cloud. The technique 700 for registration of the 3D point cloud generated by the UAV 330 can include the operator 602 providing a starting location and a heading of the UAV 330, at operation 702. The operation 702 is helpful, such as when no overhead imagery of the area or no 3D point cloud of the area is available. The heading and starting location can be determined using a compass, computing device, or the like. The heading and starting location can include an associated, estimated error. This data can be used to register the image data to the Date Recue/Date Received 2021-02-10 surface of the Earth. Using this technique, the initial heading and starting location can be used to associate the remaining 3D points to points on the surface of the Earth.

[0046] The technique 700 for registration of the 3D point cloud generated by the UAV 330 can include flying to a specified height and taking a nadir, or near nadir image of the starting location, at operation 704. The operation 704 is helpful, such as when overhead imagery or the 3D point cloud of the area are available. Overhead imagery often includes metadata, sometimes called rational polynomial coefficients (RPC), that detail locations of the pixels of the overhead imagery on the Earth. The image captured by the UAV 330 can then be registered (using normalized cross correlation of image chips, for example) to the available overhead imagery. In some embodiments, the UAV
330 can perform a LIDAR scan and take an image at the elevation. This data can be correlated with the overhead imagery to determine the starting location, at operation 706.

[0047] The geo-location registration can be performed with error (e.g., linear error or circular error, or a combination thereof) propagation. The linear error and circular error are best when correlating to a 3D point cloud, less accurate when correlating to overhead imagery, and even less accurate when only a starting heading and starting location are available. GPS data is about as accurate as an overhead imagery correlation.

[0048] The operation 704 can be helpful because the overhead imagery and the 3D point clouds are generally in nadir or near-nadir views of the tops of objects, while LIDAR from the UAV 330 sees the sides of objects in the imagery. To help correlate the LIDAR data to the overhead imagery or 3D point cloud, a nadir or near-nadir image can be generated of the starting location.

[0049] FIGS. 8-10 regard methods, systems, and devices for registering a first 3D point cloud (or a portion thereof) to a second 3D point cloud (or a portion thereof) to generate a merged 3D point cloud. One or more the first and second 3D point clouds can include an associated error. The associated error can be propagated to the merged 3D
point cloud. The error of the 3D point cloud can be used in a downstream application.
Example applications include targeting and mensuration. A targeteer (one who performs targeting) can benefit from the error to better inform their targeting location choice. A
mensuration of an object can benefit from the error as well.

[0050] The merged 3D point clouds can include error that is better than any of the 3D point clouds individually. For example, if the first 3D point cloud includes a Date Recue/Date Received 2021-02-10 lower error (relative to the second 3D point cloud) in the x and y directions and the second 3D point cloud includes a lower error (relative to the first 3D point cloud) in the z direction, the merged 3D point cloud can include error bounded by the first 3D point cloud in the x and y directions and by the second 3D point cloud in the z direction. The merged point cloud can thus inherit the better of the errors between the first and second point clouds for a specified parameter.

[0051] FIG. 8 illustrates, by way of an example, an embodiment of a system 800 for 3D point set registration and merging. The system 800 can include the processing circuitry 114 that receives tie points 808, tie point error 810, a first 3D
point set 102A, a second 3D point set 102B, and first or second point set error 812. The first or second point set error 812 includes error for at least one of the first 3D point set 102A and the second 3D point set 102B. The first or second point set error 812 can thus include error for the first 3D point set 102A, the second 3D point set 102B, or the first 3D
point set 102A and the second 3D point set 102B.

[0052] The first 3D point set 102A or the second 3D point set 102B can include a point cloud, a 3D surface, or the like. The first 3D point set 102A and the second 3D
point set 102B can include (x, y, z) data for respective geographic regions.
The geographic regions of the first 3D point set 102A and the second 3D point set 102B at least partially overlap. One or more of the first point set 102A and the second point set 102B can include intensity data. Intensity data can include one or more intensity values, such as red, green, blue, yellow, black, white, gray, infrared, thermal, or the like. One or more of the first point set 102A and the second the point set 102B can include error data. The error data is illustrated as being a separate input in FIG. 1, namely the first or second point set error 812. The error data can indicate an accuracy of the corresponding point of the point set.

[0053] The tie points 808 can associate respective points between the first 3D
point set 102A and the second 3D point set 102B. The tie points 808 can indicate a first point (xi, yi, zi) in the first 3D point set 102A, a second point (x2, yz, z2) in the second 3D point set 102B or an error associated with the tie point 808 (shown as separate input tie point error 810). The tie point error 810 can indicate how confident one is that the first and second points correspond to the same geographic location. The tie point error 810 can include nine entries indicating a covariance (variance or cross-covariance) between three variables. The three variables can be error in the respective directions (x, Date Recue/Date Received 2021-02-10 e e x xy exz y, z). A matrix representation of the tie point error 810 is provided as eyx ey eyz ezx ezy ez where the diagonal terms are respective variances in the given directions, and the off-diagonal terms are covariances between the directions.

[0054] The first or second point set error 812 can sometimes be improved, such as to be more rigorous. Sometimes, the first or second point set error 812 can be in a form that is not digestible by the bundle adjustment operation 818. The point set error 812 can be conditioned by a condition point set error operation 814 to generate an error matrix 816. The condition point set error operation 814 can include generating a covariance matrix 816 of error parameters of the first 3D point set 102A or the second 3D point set 102B. The error parameters can include seven parameters. Three of the parameters can include translation in x, y, and z, respectively. Three of the parameters can be for rotation in x, y, and z (roll, pitch, and yaw), respectively. One of the parameters can be for a scale factor between the first 3D point set 102A and the second 3D point set 102B. An example of the matrix 816 produced by the condition point set - X Xy XZ XtO X(p XK xs-yx y yz pi) ycp 3/K 3/S
ZX Zy Z ZtO Z(p ZK ZS
error operation 814 is provided as co)? coy coZ co cog) COK WS where .t is cpx cpy cpz cpco cp (plc cps KX Ky KZ KW lap K KS
_SX Sy SZ SW Sq) SK S _ translation in x, y is translation in y, where Z is translation in z, co is roll, cp is pitch, K is yaw, and s is scale.

[0055] The bundle adjustment operation 818 can receive the tie points 808, tie point error 810, first 3D point set 102A, second 3D point set 102B, and the error matrix 816 at input. The bundle adjustment operation 818 can produce a merged 3D
point set 128 and a merged 3D point set error 822 as output. The bundle adjustment operation 818 can use a least squares estimator (LSE) for registration of the first 3D
point set 102A and the second 3D point set 102B. The operation 818 is easily extendable to merging more than two 3D data sets even though the description regards only two 3D
data sets at times. The bundle adjustment operation 818 can use one or more photogrammetric techniques. The bundle adjustment operation 818 can include outlier rejection. The bundle adjustment operation 818 can determine error model parameters for the 3D data sets. Application of the error model parameters to the first 3D point set Date Recue/Date Received 2021-02-10 102A and the second 3D point set 102B, results in the relative alignment (registration) of the first 3D point set 102A and the second 3D point set 102B.

[0056] A reference number with a letter suffix is a specific instance of the general item without the letter suffix. For example, the 3D point set 102A is a specific instance of the general 3D point set 102.

[0057] FIG. 9 illustrates an example diagram of an embodiment of the relationship between ground point coordinate estimates V3 and the corresponding 3D
data set observations 12,j . In FIG. 9, three misregistered 3D data sets 902, 904, and 906 and a reference frame 924 are illustrated. First image observations 908, 910, 912 and a first associated ground point 914 and second image observations 916, 918, 920, and a second associated ground point 922 are illustrated. The ground point 914 can be determined using a least squares estimator. The least squares estimator can reduce (e.g., minimize) the discrepancy (across all observations and ground points) across all images.
The least squares estimator can project an error in one or more of the 3D data sets to an error in a registered 3D data set.

[0058] This section establishes some preliminary notational conventions and symbol definitions used for developing the formulations for the bundle adjustment operation 818. The bundle adjustment operation 818 can include identifying a ground point that reduces a discrepancy between the ground point and corresponding points in respective images, and then adjusting points in the 3D data sets in a manner that reduces the discrepancy. The term "3D data set" is sometimes referred to as an "image". For convenience, example sizes of vectors and matrices are indicated below the symbol in red. Thus, the symbol A denotes a matrix A with N rows and M columns. Column N x m vectors from R3 thus have the annotation 3 x 1. Components of a vector V are written as V = [x Y If the vector includes diacritical marks or distinguishing 3x1 embellishments, these are transferred to the components, as in = [fc 2]7' and [x' y' z]T.

[0059] Equation modeling of the relationship between points in one 3D
space to corresponding points in another 3D space is now described. A common reference space is established across all of the images. The reference space can be constructed to accommodate a simultaneous adjustment of more than two images. Correspondences Date Recue/Date Received 2021-02-10 can be formed between points in the reference space and the measured conjugate point locations in each image. The observation equation can be represented as Equation 1:
ci = (1 + s) T ( '17 ¨ i7) Equation 1 3x1 3x3 3x1 3x1

[0060] where 17 is a reference-space 3D coordinate, f/ is the observation of f7 in an image and the orientation and offset relationship between reference space and image space is taken from the orientation matrix T and offset vector V using Equation 2:
[ COCK ClASK SCOSOCK SCOSK ¨ CCOSOCK
T T (co, 0, K) = ¨COSK CCOCK ¨ SCOSOSK SCOCK CCOSOSK Equation 2 sip ¨scocip ccocip

[0061] where the symbols "c" and "s" denote trigonometric cosine and sine functions, respectively. The quantities 6' [co 0 Kr refer to rotation angles (roll, 3x1 pitch and yaw) about an image's x, y, and z axes respectively. The scalar s represents an isometric scale correction factor (nominally zero). The above form is conducive to modeling a simultaneous least squares adjustment of all images' offsets and orientations, provided that estimates of reference space coordinates for all conjugate image observations vectors are available. This form is more suitable and flexible than explicitly holding a single image as a reference for at least one of several reasons: (1) there are reference space ground coordinates that permit the potential use of ground control points, whose a priori covariances are relatively small (e.g., they carry high weighting in the solution); (2) the above formulation is suitable for a simultaneous adjustment for data that includes small or minimal overlap (mosaics), as well as, many images collected over the same area (stares) or any combination in between;
and (3) a single image can effectively (e.g., implicitly) be held as a reference by appropriate a priori weighting of its error model parameters.

[0062] The symbol 17 will be referred to as a ground point (akin to tie point ground locations and ground control point locations in a classical photogrammetric image adjustment). The symbol V will be referred to as a ground point observation (akin to image tie point observation locations in a classical photogrammetric image adjustment).

[0063] Unlike the classical photogrammetric treatment, 17 and -17 are both "on the ground" in the sense that they both represent ground coordinates in 3D (in the classical imagery case, the observations are in image space and are thus 2D
coordinates). Further, the point may very well not be "on the ground" but could be on a building rooftop, Date Recue/Date Received 2021-02-10 treetop canopy, etc. However, the terminology "ground point" and "ground point observation" will be used.

[0064] If] is taken to be the index of an arbitrary ground point and i to be the index of an arbitrary image, the observation equation (Equation 1) can be written as Equation 3 = (1 + si)TA ¨ Vi) Equation 3

[0065] where Vi rii jij j1 as the jth ground point, Vi [xi yi zilT as the offset vector between image i and the reference space origin, and where COiCKi C(DiSKi StOiSOiCKi S (DiSKi ¨ C(DiSOiCKi Ti T GO i, 0i, Ki) = ¨CO iSKi C(DiCKi ¨ SO) iSOiSKi S (DiCKi C(DiSOiSKi [
SOi ¨SCOiC0i CCOiCO i Equation 4

[0066] is the orientation matrix between image i and the reference space frame and where si is image i scale correction factor. Thus, Vii is the coordinate of the ith image's observation of ground point].

[0067] If a particular ground point is found in two or more images, it can serve as a point which ties the images together (one of the tie points 808). These are generically referred to as tie points (or tie points). A single tie point is often referred to as a collection of image observations (with coordinates) of the same point on the ground along with the corresponding ground point (with coordinates).

[0068] Since the observations over many images i are the measurements containing error, the true ground point coordinates are generally unknown. To facilitate this, an initial estimate of the ground point location can be computed. The initial estimate is provided as Equation 5 Ei -1. Equation 5 1 11_11 El j ij

[0069] The ground points themselves are treated as derived (but unconstrained) observations and allowed to adjust in performance of the operation 818. There can be an observation of interest whose true ground coordinates are well known. These are classically called ground control points (or GCPs). Since this development can accommodate both GCPs and tie points, the more general terms of "ground point"
and "ground point observation" are sometimes used (as contrasted with "tie point ground coordinate" and "tie point observation").
Date Recue/Date Received 2021-02-10

[0070] The bundle adjustment operation 818 can operate on two or more images taken over a same area (with observations for tie points, sometimes called a stare scenario); two or more images taken in strips (forming a mosaic of data, with 2-way, 3-way, or m-way observations in strip overlap regions); tie points in which the corresponding ground points may appear in two or more images, incorporation of GCPs for features in imagery, providing an absolute registration; accommodation of a full covariance for tie point observations. This is conducive for tie point correlation techniques which are highly asymmetrical (e.g., as long as the asymmetry can be characterized as a measurement covariance).

[0071] The relationship between ground point coordinate estimates 12] and the corresponding image observations Vii can be understood as a stare scenario between three misregistered images.

[0072] For the development of the LSE formulation (and associated preprocessing) that can be performed by the bundle adjustment operation 118, more definitions are provided in Table 1.
Table 1: Definitions of Symbols Symbol Definition VRw Location of the origin of the reference frame with respect 3x1 to the world frame. This is thus the location of the reference frame coordinatized in the world-frame.
Translation of ith image with respect to reference frame 3x1 origin.
17i ei Orientation angles of it" image with respect to image 3x1 frame origin 01 Oi KIlT
si Isometric scale factor correction for the ith image.
tx oco Initial (zeroth-iteration) value for i7, thus assumed to be 3x1 with respect to reference frame origin.
c') Initial (zeroth-iteration) value for Oi. Each element is 3x1 taken to be zero.
si( ) Initial (zeroth iteration) value for si. Nominally si(0) 0.
tx Date Recue/Date Received 2021-02-10 Symbol Definition Orientation matrix for ith image built from Oi 3x3 Ground point coordinates for ground point j with respect 3x1 to the reference frame origin [j2i 2l1T
0.0) Initial (zeroth-iteration) estimated value for 12./
v 3x1 Ground point observation coordinate of ground point j on 3x1 image i coordinatized in the world frame (e.g., these are UTM coordinates of the ground point observation location).
Ground point observation coordinate of ground point j on 3x1 image i. These are implicitly assumed to be coordinatized in the local image frame for image i.
2illT
A A priori covariance of it h image translation, orientation 7x7 and scale correction parameter vector 147i A priori parameter weight matrix for image i.
7x7 =i-'WI
A priori covariance of ground point j 3x3 1;17- A priori weight matrix for ground point j.
3x3 147i =
A priori covariance for observation of ground point j upon 3x3 image i j A priori observation weight matrix for observation of 3x3 ground point] upon image i. W11 = 1 General Indexing Number of images Image index. i E {1,2, ===, Total number of ground points Ground point index j E 02, ===, Date Recue/Date Received 2021-02-10 Symbol Definition Total number of ground point observations Ground point observation index b E 02, === , (19) Non-linear least squares iteration index Gi The index set of all ground points appearing in image i.
Thus Gi {1,2, === ,n}
O The index set of observations of ground point] over all images.
Thus Of g f1,2, = = = , q}
The index set of images upon which ground point] is an observation. Thus g ti,2,===,m}
Mg Mapping of observation index to Ground point index.
mg gives the ground point index (E t1,2, = = = , np for a specified observation index b E f1,2, = = = , q}.
151 Cardinality of set S (e.g., the number of index elements in set S).

[0073] Ground point observations can be indexed by ground point] and image i (as in PO or by linear indexing, b (as in Vb). Use of the subscripting depends upon the context. In the former, it is of interest to characterize the fact that a particular ground point] appears on a particular image i. In the latter, it is of interest to enumerate all observations independent of to which image or to which ground point they refer.

[0074] Since some 3D point set data is presented in a "world" space coordinate system (e.g., Universal Transverse Mercator (UTM) map projection) and since the observation Equation 3 is image dependent, some coordinate frame definitions and transformations can aid understanding.

[0075] If it is assumed that ground point observation locations are specified in world coordinates, it is of interest to transform the ground point observation locations to be "image" relative. Further, it can be of interest to obtain the ground locations and image offsets themselves to be relative to a "local" reference coordinate frame.

[0076] A motivation for a local reference coordinate frame can be to remove large values from the coordinates. For example, UTM coordinates can typically be in the hundreds of thousands of meters. This makes interpretation of the coordinates more difficult, for example, when examining a report of updated coordinate locations. A

Date Recue/Date Received 2021-02-10 motivation for an image-relative coordinate frame can be so that the interpretation of the orientation angles comprising the Ti matrices can be relative to the center of the data set.
This is contrasted with the origin of rotation being far removed from the data set (e.g., coincident with the local reference frame origin in the mosaic scenario).

[0077] In both cases, the transformations between coordinate frames simply involve a 3D translation. The mnemonics W, R and I are used to denote the "world", "reference" and "image" coordinate frames, respectively. To facilitate the transformations, the following convention is established. A superscript on a vector denotes the coordinate frame to which it is referred. Thus 7Y corresponds to the world space coordinates of a particular tie point observation, while I7i).1 and represent the same tie point observation but referred to the reference frame and image frame, respectively.

[0078] Following the above convention, the symbol can VI represent "the location of the origin of frame A coordinatized in frame B". Thus, VRw can represent the location of the reference frame in the world coordinate system (e.g., UTM
coordinates of the origin of the reference frame). The relationship between an arbitrary vector VR
coordinatized in the reference frame and the same vector Vw coordinatized in the world frame can be represented by Equation 6 Vw = VR + 1/Rw Equation 6

[0079] The reference frame can be established as an average of all of the world-space coordinates of tie points. This offset (denoted Vi') can be determined using Equation 7

[0080] Vw --iEq R b=1 b Equation 7

[0081] For simplicity, it can be assumed that the reference frame origin, referred to by the world frame, can be computed by a process external to the bundle adjustment operation 818 (e.g., by the process that assembles the tie points 808 for use in the bundle adjustment operation 818).

[0082] The image frame (e.g., a frame defined on a per-image basis) can be the world coordinates of the center of an image. Under the assumption that there are bounding coordinates in the image data (specifying the min and max extents of the data in world-frame X, Y and Z), the center of the data can thus be taken to be the respective averages of the min and max extents. Since this image frame refers to world space, the computed offset is denoted Viw. If bounding coordinates are not available, value for Viw Date Recue/Date Received 2021-02-10 is taken as the average of the tie point locations over the specific image i, as described in Equation 8 Viw = E JEG tj Equation 8 IGti t

[0083] The image frame offset in reference space coordinates is taken to be the initial value for V(0) on a per image basis. Thus, for each image i, an external process can compute reference frame coordinates according to Equation 9 = 171w ¨ VRw Equation 9

[0084] Since the tie point observation values can be input in world coordinates and since the observation equation domain assumes reference frame coordinates, some preprocessing of the input data can help make it consistent with that assumed by the observation equation (Equations 1 or 3). The tie point observation coordinates can be converted from world space to reference space. This can be performed for each observation per Equation 10.
= Pity ¨ VRw Equation 10

[0085] Next, since the true ground point coordinates used in Equation 3 can be unknown, they can be estimated. The ground point coordinates can be assumed to be coordinatized in the reference frame. The initial estimated values for the ground coordinates of each tie point can be computed as an average of the ground point observations over all images in which it appears as described by Equation 11 Equation 11 1 11_11 El tj

[0086] Since the true locations of the tie point ground coordinates can be treated as unknown, the a priori covariance can reflect this by treating the errors in the ground coordinates as numerically unconstrained (in units of meters squared) as described by Equation 12 diag([1012 1012 10121) Equation 12

[0087] The tie point observation coordinates for use in the observation equation can be converted to image-relative coordinates using Equation 13.
= ¨ Vi¶)) Equation 13

[0088] Next, a least squares formulation and solution are discussed.
Since the observation equation, Equation 1 or 3, is non-linear in the orientation angles that form Ti, the least squares problem becomes a non-linear least squares problem.
Equation 3 Date Recue/Date Received 2021-02-10 can be linearized. Solving the linearized equation can be a multidimensional root finding problem (in which the root is the vector of solution parameters).

[0089] For simplification in notation of the linearization, consider a fixed image and a fixed ground point. Let the unknown error model parameters (offset, orientation, and scale correction) be represented by Equation 14:
X [x y z cp KS]T Equation 14 7x1

[0090] The observation equation for ground point observation V can be written as Equation 15 = (1 + s)T(I7 ¨17) Equation 15

[0091] where T is the true image orientation matrix, I/ is the true image translation vector, 17 is the true ground point coordinate and -V is the corresponding ground point observation coordinate.

[0092] If one wishes to include the ground point coordinates 12 as additional observations, the solution for X and 17 can be cast as a root solving problem based on Equation 16 F(X; V) = 0 Equation 16

[0093] where F(X; V) = ¨ (1 + s)T(V ¨17) Equation 17

[0094] In vector font', the function, F, can be represented by Equation b2I = ¨ (1 + s)T ([91 _ ki) Equation 18

[0095] The function F can be approximated using a first-order Taylor series expansion of F about initial estimates X(') and V( ) as in Equation 19 a F(0) a F(0) 7 0 = F (X ; 17) F (X( ); i7()) + ¨ A A Equation 19 a x a v

[0096] where X((0 is an initial approximation of X, j7(0) is an initial a F(0) a F(0) approximation of 17 , the Jacobians ¨ax and ¨ay are the partial derivatives of F
evaluated at X( ) and V(0) respectively, A [M Ay .62 Act) Alp AK AS]T is a vector of corrections to X, Aji .62] is a vector of corrections to V. The values for X((0 and 17(0) are discussed in Table and in sections Error!
Reference source not found. and Error! Reference source not found..

[0097] The Jacobians can be written as Equations 20 and 21 Date Recue/Date Received 2021-02-10 -af1- -af1 afi af1 afi af1 afi afi-=axa ay .92 .90) av ak as = aF af2 af2 af2 af2 af2 af2 af2 af2 B ¨¨Equation 20 3x7 ax ax ay, .92 .90) av ak as a I's afs a I's a I's a I's afs a I's a I's -ax- a ay .92 .90) av ak as _ -a ft- -=all all afi-a Va a9 a2 aF af2 af2 af2 af2 ¨= ¨= ¨ -- --Equation 21 3x3 av av a9 a2 af3 afs a I's afs _aR a9

[0098] Note that the dot symbols are merely notations, following the classical photogrammetric equivalent, and do not intrinsically indicate "rates," as is sometimes denoted in other classical physics contexts.

[0099] In matrix notation, Equation 19 can be written as pm] [Al = _F(x(0); i7(0)) Equation 22 ax av 2j

[00100] or 3x7 3x3 7A
[ I[ Ax1 E
_ 3x1 Equation 23 3x1

[00101] Since the problem is nonlinear, the estimation of the parameter vector can be iterated (via a multi-dimensional extension of the Newton-Raphson method for root finding, or other technique). The solution can include re-linearization at each iteration. The re-linearization can be performed about the most recent estimate of the parameter vector. The linearized form of Equation 22 at iteration (p) can be represented as in Equation 24.
,N a F(p) a F(p) 7 F(X; V) F (P) ; (Pi +x Equation 24 av

[00102] where X(P) is the pth iteration estimate of the parameter vector X, (P) is ap(p) the pth iteration estimate of 17 , ¨a x is the Jacobian of F with respect to X
evaluated at is the Jacobian of F with respect to 12 evaluated at 12(P), A is a vector of corrections to X for the p' iteration, and 21 is a vector of corrections to 17 for the p' iteration.

[00103] With each iteration, the parameter and ground point vectors can be updated with the most recent correction as in Equations 25 and 26.
X(P) = + A Equation 25 Date Recue/Date Received 2021-02-10 17(P) = 17(3-1) + A .. Equation 26

[00104] .. For the initial iteration, initial values for X(') and I7( ) can be estimated as discussed previously. The system represented by Equation 24 is now linear in A and A. A linear solver can be used to solve for the parameters.

[00105] .. For a particular image i and a particular ground point j, Equation 23 can be written as Equation 27 Ai I
[ [E11 Be.] 7x1 ,..

3 Equation 27 1-3x7 3x3-1 [14 x ¨ ;--)1<j

[00106] .. The discrepancy vector for the loth iteration is thus be represented as in Equation 28 Ei J)) ¨ ¨ F(X(P); I7(P)) .. Equation 28 ¨

[00107] .. and thus (p) iv _ = ¨1_ ii (1 + si)TA ¨17i)] Equation 29

[00108] .. To accommodate a simultaneous solution of all images and ground points, Equation 27 can be extended as - A,-- En o ¨ o PH o ¨ o - A2 En B12 0 === 0 0 ij12 ¨ 0 i E12 0 =.0 A. i = Equation 30 0 B21 ¨ 0 :621 0 i 0 A1 E21 _ 0 0 === Bmõ 0 0 .. === Bmn--An -

[00109] .. Equation 30 can be re-written as Equation 31 BA = E Equation 31

[00110] .. then the normal equation matrix can be represented as Equation 32 or Equation 33 (BTB)A = BTE Equation 32 ZA = H Equation 33

[00111] .. It can be less efficient to form B as in Equation 30, for one or more of the following reasons: (1)B is very sparse; (2) the quantities Eii and Bii are nonzero if and only if ground point] is observed on image i. For this reason, the classical development of the normal matrix BTB and right-hand side vector BTE uses summations over the Date Recue/Date Received 2021-02-10 appropriate indexing. These summations are provided in the normal matrix partitioning below.

[00112] The foregoing equations form a foundation for the present problem that is sufficient for development of the normal equations, examination of the normal matrix structure and formulation of the normal equation solution.

[00113] The normal equation can be written as in Equation 34 ZA = H Equation 34

[00114] The matrices can be partitioned as in Equations 35-37 -F Vi/ i -I
Z = _ Equation 35 I_ NT + WI
A =[A-1 A Equation 36 H = [Wel Equation 37 IKI WC

[00115] The quantities k,k, e and C are described in more details elsewhere herein.

[00116] Combining Equations 35, 36 and 37 yields Equation 38 NT +N V171[11 = [kij [WW C1 Equation 38

[00117] The matrix Z can thus be represented as Equation 39 Z = N +W = 7mx7m 7mx3n 7mx7m Equation 39 3nx6m 3nx3n 3nx3n

[00118] The matrix /V can be written as Equation 40 7x7 N = 7x7 Equation 40 7mx7m 0 0 Nn, 7x7-

[00119] and analogously W can be written as Equation 41 7x7 W = 7x7 Equation 41 7mx7m 0 0 Wm 7x7-

[00120] The block entries of Ni can be defined as in Equation 42 Date Recue/Date Received 2021-02-10 N i =EjEGiETjWijEij Equation 42 7x7

[00121] The subscripts ij on the Pi./ matrices indicate that they are a function of image i and ground point j.

[00122] The matrix N can be expanded as in Equation 43 0 0 -3x3 = 3x3 Equation 43 3nx3n 0 0 Nn 3x3-

[00123] W can be expanded as in Equation 44:
0 0 -3x3 W = 3x3 Equation 44 3nx3n 0 0 Wn 3x3-

[00124] The block entries of Equation 43 can be defined as in Equation 45 N1 = EiEo; BWjBj Equation 45 3x3

[00125] The matrix N from Equation 39 can be expanded as in Equation 46 NIA N12 = = = Nln 7x3 7x3 7x3 N21 N22 = = = N2n N 7x3 7x3 7x3 Equation 46 7mx3n Nmi Rm2 = = = kin - 7x3 7x3 7x3 -

[00126] The block entries of N from Equation 45 can be defined as in Equation Nii = Pri Wu Dui E 1, = = = , m), j E 1, = = = , n) Equation 47 7x3

[00127] In a similar development the right hand side matrix H from Equation 34 can be expanded as in Equation 48 H= kõ? n _[W [e 1 Equation 48 K1 0 Wild Date Recue/Date Received 2021-02-10

[00128] The subblocks of H can be defined as in Equations 49 and 50 ki = EJEGtfiritil7 i;Ei; Equation 49 7x1 kj =TitiViiEii Equation 50 3x1

[00129] with the discrepancy vector Eii defined as in Equation 29.

[00130] e(P) = ¨ e(0) Equation 51

[00131] C(P) = - C(13) Equation 52

[00132] The values for e(0) and ;( ) are the initial parameter values. The initial values for the translation parameters portion of e(0) can be taken to be the V() as computed in Equation 9. The initial values for the rotation parameters portion of e(0) can be taken to be zero.

[00133] The initial values of CO) can be taken to be the values of the ground point coordinates 17./( ) as computed in accord with Equation 11.

[00134] The solution to the normal equation matrix on iteration (p) can be determined as in Equation 53 .6(P) = Z-111 Equation 53

[00135] At each iteration, the parameters can be updated via Equations 51 and Equation 52 and the normal matrix can be formed and solved again. The process can continue until the solution converges. Examples of the convergence criterion can be discussed in the following section.

[00136] Since the solution is iterated, a convergence criterion can be established.
An example of a convergence criterion is to compute the root-mean-square (RMS) of the residuals as in Equation 54 R(P) = ss Equation 54 q-7m

[00137] The value in the denominator of Equation 54 represents the number of degrees of freedom (e.g., the number of observation equations minus the number of estimated parameters).

[00138] Since typically q >> 7m Equation 54 can be estimated as in Equation 55 R(P) = Equation 55

[00139] The condition q >> 7m can be guaranteed with sufficient redundancy of ground point observations as compared with the number of images (e.g., enough tie Date Recue/Date Received 2021-02-10 points are measured between the images so that the aforementioned condition is satisfied).

[00140] Convergence happens when the residuals settle to the same values on consecutive iterations. The convergence criterion can be IR(P) R(P-1-)1 <8 Equation 56

[00141] where 6 is a prescribed tolerance.

[00142] A rigorous formulation for the standard error of unit weight (to be used in error propagation discussed elsewhere) is provided in Equation 57 2 En¨iEtE/ W E Ein eTw +El' eTiv jCj [0-03)] _ ndof Equation 57

[00143] where ndof is the number of degrees of freedom¨the number of observation equations minus the number of error model solution parameters:
ndof = q ¨ 7m Equation 58

[00144] Since blundered points can be effectively removed from the solution via deweighting, the number of observations remaining effectively doesn't include the blunders. To be strictly correct, the value for q in Equation 58 can be the number of non-blundered observations.

[00145] The full form of the matrix Equation 34 can be reduced under the assumption that the errors in the ground point locations are uncorrelated.
Under this assumption, the error covariance matrix of the ground point locations becomes a block-diagonal matrix of 3 x 3matrix blocks. Since it is a sparse matrix, its inverse is easily computed by inverting the 3 x 3 diagonal blocks. The development in this section reformulates the normal equations taking advantage of this. The result is a reduced normal equation matrix in which the size of the normal matrix is 6m x 6m instead of (6m + 3n) x (6m + 3n). This gives the obvious advantage that the size of the normal matrix is much smaller and remains invariant with the number of ground points.

[00146] The reduced system formation is sometimes referred to as a "ground point folding," since the ground point portion of the reduced normal matrix is incorporated into the image portion. The development of the reduced normal equation begins with the original normal equation from Equation 34 and repeated as Equation 59 ZA = H Equation 59

[00147] To facilitate ground point folding into a reduced normal equation matrix, Equation 59 can be re-written as Equation 60 Date Recue/Date Received 2021-02-10 [)2<7m 72<3n1 [7m A x11 = [71i: 7m I
x11 Equation 60 3nx6m 3nx3n 3nx1 3nx1

[00148] where 2 = N + Vi7 Equation 61 2 = N Equation 62 2 = N + W Equation 63 H = k ¨ we Equation 64 ii = k - V17' Equation 65

[00149] Suppose a matrix system ZA = H is partitioned into blocks of the appropriate sizes as [C DABI [ c/ abl = [c1 Equation 66

[00150] where the matrices A and D are both square.

[00151] Further, assume that matrix D is non-singular and can be represented as a sparse block diagonal matrix. Then [A ¨ BD-1-C] [a] = [c ¨ BD-1d] Equation 67

[00152] Applying Equation 67 to Equation 59 provides the reduced normal matrix equation [2 - 22-1-27][A] = [fl - 22-111] Equation 68

[00153] The reduced normal equation matrix can be written as in Equation MA = C Equation 69

[00154] where M [2 ¨ 22-121 and C [1-*/ ¨ 2 2-11-1].

[00155] Next it is of interest to examine the form of the components of the reduced system for an efficient implementation. Let 2 22-1-2T . Then -21-1 0 === 0 -(7x3)... (<731) (3)<03) 2-1 2r1 ... 427'1 === 0 (3x7) (3x7) 2 = i ... i (3x3) Equation 70 7mx7m -Zmi === 2,2n-'in = = = 2Tnii, (7x3) (7x3) 0 0 === 27T1 (3x7) (3x7) (3x3)-

[00156] By extension - ZY=1 2 ' 1f2j-12f Z Y=1 2 -1 ' if Y
... Z=1 2 -1 T
if 2f2172 i 2 = ZY=1 22i2i 14 zY=1 22j2j 121f === ZY=1 22i 2-i 127nTi Equation 7mx7m -rf1=1 2172/2f-12rf ZY= 1 2172/2f-12'f ... ZY=1 2m/2f-12172T i Date Recue/Date Received 2021-02-10

[00157] The blocks of 2 in Equation 71 can be the equivalent Nii as defined in equation 47.

[00158] The assumption that errors in the a priori ground points are uncorrelated yields Equation 72 1,71 0 0 -3x3 0 1-1 === 0 W =o 32x3 Equation 72 3nx3n 0 0 === n-1 3x3-100159] where IT' is the inverse of the a priori covariance matrix for ground point]. Thus 27' = [N. + Equation 73 3x3 [00160] The general row and column term for 2 can then be given by 2r,c = Erj1=12rj2j 12 cT j Equation 74 7x7 [00161] and, by the definition of Z1, 2r,c is zero if and only if images r and c have no ground points in common. Also note that 2 is block symmetric. Thus, in its formation, only the upper block triangle need be formed, followed by reflection of the upper right triangle to the lower left triangle for completion of the matrix formation.
[00162] The matrix M can thus be written as in Equation 75 - 21 0 === 0 -(7x7) 0 z2 0 M = (7x7) (7m x 7m) = =
0 0 ===
(7x7)--E7=121;2;-12; =1- 1] Ey=121;2;-12,,T
vr!. vr!. vr!.
L-,j=1'=2j''j L-,J =1 j ]=1 j mj Equation 75 vr!. 2 .2712T
-L-,J =1 '-' ITU j'-' 1] L-,J =1 '-' ITU j'-' 2] =1 In1 j 1n1 [00163] The reduced matrix M can be formed by first storing the diagonal entries of 2 and then subtracting the summed entries of the subtrahend in Equation 75 (namely the 2rx defined in Equation 74).
[00164] Since the subblocks of the subtrahend are merely summations over the ground point indexes, j, the matrix, M, can be built by iterating over the ground points Date Recue/Date Received 2021-02-10 (assuming the minuend of Equation 75 on-diagonals were formed in advance) and subtracting out the contributions for a particular ground point in the appropriate place within M.
[00165] The constant column vector C can be formed similarly with some of the same matrices:
c, . [ 1 ri,=1 : _ 11 Eill 2 .27111.
1 In1 1 1 Equation 76 [00166] After the matrices M and C are built, the solution vector for the adjustable parameters from the reduced system can be computed as A = M-1C Equation 77 [00167] The solution vector can be decomposed into per-image-adjustable vectors Ai for each image i as in Equation 78:
A, A ,Zi.)21 Equation 78 A
[00168] After the solution vector A for the image-adjustable parameters is obtained, the solution vector A for corrections to the ground point positions can be extracted (or -unfolded") from the reduced system. To formulate the extraction, Equation 67 can be used to obtain Equation 79 A = 2-1[Ii ¨ 2T A] Equation 79 3nx1 [00169] If -Ai -3x1 A = 3x1 Equation 80 3nx1 An -3x1-[00170] represents the correction vector for the ground points then -A1- - 21-1 o === 0 - - a - -3x1 (3x3) 3x1 2T1 === 2inT1 7x1 A2 0 22i- === 0 112 (3x7) (3x7) A2 ..
3x1 = (3x3) 3x1 ¨ : . 7x1 Equation : =. : in ¨ 21Tna An 0 0 === 2n-1- tin (3x7) (3x7) Ain -3x1- _ (3x3)- _-3x1- -7x1-Date Recue/Date Received 2021-02-10 [00171] where A iis the adjustable parameter correction vector for image i. Thus A. = 271[1:1j ¨ EiEli TiAi] Equation 82 3 x 1 [00172] where I is as defined in Error! Reference source not found. (the index set of images upon which ground point j is an observation).
[00173] This section provides formulations for extraction of a posteriori error covariances for ground points. If a priori sensor model error estimates are available (and reliable), the errors may be propagated to the space of the registration error models. In this case, the error propagation is a rigorous predicted error for the accuracy of the a posteriori ground point locations.
[00174] The a posteriori error covariances of the image parameters are the appropriate subblocks of the inverse of the reduced normal matrix M-1 from Equation 69 (after application of the variance of unit weight, as described at the end of this section). For the full normal matrix solution, the a posteriori error covariance can be the inverse of the normal matrix, Z-1, times the variance of unit weight. For the reduced system, however, the a posteriori error covariances of the ground points can be extracted from M-1 by unfolding. To facilitate this, the full normal matrix can be written as (7m x 7m) (7mx 3n) Z = Equation 83 (3n2xT7m) (3n2x 3n) [00175] Denote the inverse matrix blocks as Z-1 [ Equation 84 IT I
[00176] Note that, and I as defined are distinctly different from those defined in previous sections. (The symbols in the present section are a posteriori covariances and those in previous sections are a priori covariances). However, this subtle distinction is not problematic if the appropriate context is adhered.
= + 2-12T 22-1-Equation 85 3n x 3n [00177] The a posteriori covariance between ground points r and c can be represented as block element rx oft. With n as the number of ground points and m as 3x3 the number of images, Date Recue/Date Received 2021-02-10 Ill = = = Iln [
=
Inl = = = Inn 21-1- 0 ==== 0 0 2-'=== 0 0 0 === 277,1-2V- 0 === 0 21- 0 === 0 0 2ii- = = = 0 VII 0 221- === 0 0 0 ¨ Z.2w-'¨ 27T1- Equation 86 [00178] The rth row of 1 involves only 47' of the first2-1 matrix in term two of Equation 86. Similarly, the cth column of 1 involves only 2,71 of the second 2-1- matrix in term two. Thus frx = 8(r, c)2,71- + 4712T l'22,71- Equation 87 [00179] where the delta function can be f1 if r = c 8(r, c) Equation 88 (0 otherwise [00180] Now the form of the (r, c) block of 27'11 is derived. Let E1 1'12 = = = Elm (7x7) (7x7) (7x7) 211 === 21n (3 (37) /21 /22 === I2n2 (7x3) (7x3) G 27'11 A = i '.. : (7x7) (7x7) (7x7) 3inx3in 2in ¨ 2Tnris i i *** i Ll ¨
2inli, (3x7) (3x7) An1 Im2 === LT (7x3) (7x3) -(7x7) (7x7) (7x7)-Equation 89 [00181] The rth row of G involves only the rth row of 27' and the cth column of G involves only the cth column of Z. Thus G r,c = [ZT411 = = = Zinr2mi Zrõ.212 + = = = + Zmr2m2 === 1Z
3x3 Zic = [Esm-1 nr2s1 Esm-1 nr2,32 [ === 1 Z`' Equation 90 zinc m m = 1 [ 1 nr 2st] Ztc t=1 s=1 [00182] Now 2ii = 0 if ground point j is not an observation on image i.
[00183] Thus Tr,c = EtedEsElr 2 TrEst] 2 tc. Equation 91 [00184] where b is the index set of images upon which ground point j is an observation. Substituting Equation 91 into Equation 87 yields Date Recue/Date Received 2021-02-10 1r ,c = 8(r, c)2 c 1 2 rl[EtElc[EsElr 2 sTrkti2tc12 c 1 Equation 92 [00185] The a posteriori covariance is usually defined by scaling the inverse of the normal matrix by an estimate of the variance of unit weight. An estimate of the variance of unit weight is denoted as [o-(P)]2 and is provided in Equation57.
Thus, the above formulation can be used, but instead defining E 1 [o-(P)]2 [Z1 ¨ Equation 93 LT 11 ¨ ¨
[00186] For a full normal matrix solution, Z-I- is readily available, thus the a posteriori covariance of the error model parameters and ground points can be the right hand side of Equation 93.
[00187] The right hand summand of Equation 92 includes the factor WO]2 since it includes Est. However, the left hand summand does not include the factor.
This can be compensated for by a modified form of Equation 92 1r ,c = [a(P)]28(r, c)2c 1 + 2r1[EtElc[EsElr 2 sTrE st12 tc12 c 1 Equation [00188] If the standard error of unit weight o-(P) is deemed to be unreliable (e.g., is much greater than unity) this may be an indicator of improper (or incorrect) a priori error covariance in the process. One can still, however, be able to provide a reliable error estimate from the least squares process by simply forcing the standard error to one (e.g., by setting cr(P) <¨ 1.0 in Equations 93 and 94.
[00189] FIG. 10 illustrates an example of an embodiment of the operation 818.
The operation 818 can include 3D data set registration with error propagation.
The operation 818, as illustrated, includes initializing solution and corrections, at operation 1002; determining discrepancies, at operation 1004; determining the normal equation, at operation 1006; updating parameters based on the determined normal equation, at operation 1008; determining discrepancies, at operation 1010; determining error, at operation 1012; and compensating misregistration of the first 3D point set 102A and the second 3D point set 102B, at operation 1014.
[00190] The operation 1002 can include setting the solution vector X and the correction vector AX to the zero vectorl:
X = 0 7x1 AX= 0 7x1 Date Recue/Date Received 2021-02-10 [00191] The solution vector X can be set to a fixed-point location for the linearization. If an a priori estimate is available, it can be used here in place of the zero vector.
[00192] The operation 1004 can include computing the discrepancy vector for each observation as provided in Equation 29. The operation 1006 can include building the normal equations matrices and solving for the correction vector as provided in Equation 53. The operation 1008 can include updating the parameter vector for the current iteration as provided in Equations 51 and 52. Details of the operation 1008 for unfolding of the ground points for the folded normal equation solution is provided via pseudocode below.
[00193] The operation 1010 can include computing the discrepancies as provided in Equation 29. The operation 1012 can include computing a standard error of unit weight via Equation 57. Note that the standard error can the square root of the left hand side of the Equation 57 (e.g., o-(P) = -I[o-(P)]2).
[00194] If the delta between the current and previous standard error of unit weight is less than the convergence criterion in absolute value, the solution nominally converged. To accommodate blunder rejection, the convergence criterion check can be augmented with a check to see if the blunder weights should be used in continuation of the solution (-useBW", indicating to use "blunder-checking weighting"). If convergence occurs and useBW is true, this is an indicator to perform blunder checking, and this time using a normalized residual computation in order to check for blunders on the next iteration.
[00195] If useBW is true, blunders can be computed. If there are blunders remaining, the blunder "cycle" number is incremented and the process is repeated with the correction vector reset to a priori values (e.g., go to operation 1002).
If there are no blunders remaining, a check can be performed to see if the number of post convergence blunder cycles can be set to zero. This check can be performed to effectively force one more solution after all blunders have been eliminated.
[00196] If useBW is false and it is currently a first iteration of a blundering solution, useBW can be set to true. This has the effect of forcing the normalized residual blunder iteration for determining the blunders on subsequent iterations. In this case, a solution has converged but normalized blunder residuals have not been computed.
Setting useBW to true can forces this to happen on the next solution iteration. The Date Recue/Date Received 2021-02-10 solution can be iterated by going to the operation 1006. If there are no more blunders and the number of blunders is not zero, this indicates the "non-blunder iteration"
solution has converged.
[00197] The operation 818 can include providing a report that includes an iteration number, current correction vector AX, current iteration estimates of parameters and ground points (e.g., as computed in equations 51 and 52), standard error of unit weight (e.g., as provided in Equation 55). The operation 818 can include a check for non-convergence by examining the current iteration number with a maximum number of iterations, M. If the number of iterations exceeds the maximum, stop the iteration process. The solution did not converge. An exception can be raised and the operation 818 can be complete.
[00198] The following is a pseudocode outline for the operation 818 for computing the full normal equations solution. This first pseudocode does not include ground point folding.
[00199] The names in {braces} ce s } allude to method (e.g., function) names in a software implementation. Also, within the development below, ground point indexes, ground point observation indexes and image indexes are assumed to be zero-relative.
For efficiency in the implementation, the following elements can be cached in a per-image workspace object, which is updated with each iteration:
1. Values for trig functions sco, cco, sO, CO, SK, CK as given in Equation 4.
2. Three 3 x 3 partial derivative matrices of T with respect to the three angles as given in the following Equations 95-97:
0 ¨StOSK CGOSOCK CGOSK StOSOCK
aT
¨ = 0 ¨StOCK ¨ CGOSOSK CtOCK ¨ StOSOSK Equation 95 .90) 0 ¨mop ¨scoop ¨SOCK stocOcK ¨ctocOcic TO ¨aT = sOsic ¨scocOsic ccocOsic Equation 96 acp scos0 ¨ccosip aT
¨Cc/SK CtOCK ¨ StOSOSK StOCK CGOSOSK
T' = ¨COCK ¨CCOSK ¨ StOSOCK ¨StOSK CGOSOCK Equation 97 3. Rotation matrix Ti as given in Equation 4 4. 3 x 1 vector 17i [00200] The pseudocode begins by setting the non-linear least squares iteration index (p) to zero.
Date Recue/Date Received 2021-02-10 1. {initializeData} Store nominal initial values for the solved parameters for each image i E tO,===,m ¨ 1). The initial values for Vi( ) and 0(0).
[00) Equation 98 s.(0) 71x1 2. {initializeData} Compute a. Initial ground point coordinates via Equation 11 for each ground point j E
tO,===,n ¨ 1). These form C./M.
b. Initial image cache data as described previously.
3. {outputIterationReport} Output initial iteration report (input parameters, initial ground point and ground point observation coordinates.) 4. {initializeMatrices} Block partition the R portion of Z into m subblocks, each of size 6 x 6.
5. {initializeMatrices} Block partition the N portion of Z into n subblocks, each of size 3 x 3.
6. {initializeMatrices} Block partition H similarly 7. {computeAndStoreDiscrepancies} For each ground point index j E
to, ===, n ¨ 1) and each observation of that ground point i E Op compute the discrepancy vector Eij given in Equation 29 as follows.
for] E f0,===,n ¨ 1}
{
Fetch the most recent ground point position 17j.
for i E b and observation index b E Of {
a) Retrieve image-cached values for Ti and Vi b) Retrieve ground point observation Vij = Vb c) Apply the observation equation to obtain the projected value for 17j.
= (1 + si)TA ¨ Vi) Equation 99 d) Compute and store the discrepancy vector for observation b as in Equation EL, Eii = ¨[Vii ¨ Vii] Equation 100 end for i Date Recue/Date Received 2021-02-10 Wend for j 8. Compute the standard error of unit weight as in Equation 57.
9. {buildNormalMatrices + initialize Weights} Zero Z and H and initialize Z (and likewise H) with the block weight matrices on the diagonal.
This involves setting the blocks of Z to the subblocks of W and W. and setting the subblocks (subrows) of H to ¨we and ¨WC.
10. {sumInPartials} Loop over ground points and images containing the ground points and sum in the contributions of the .6' and .6 matrices into Z and H.
for] E 0, = = = , n ¨
for i E
Retrieve Eb Eii as computed in Equation 100 Compute Bii and B11 as in Equations 101 and 102.
[(1 + si)Ti (1 + si)Aii ¨TiYiii Equation 101 3x3 3x3 3x1 3x7 = ¨(1 Si)[Ti] Equation 102 3x3 3x3 Retrieve observation weight matrix Wi : Sum into Z.block(i, i) Nj : Sum Pi; into Z.block(m + j, m +j) Nii : Sum PiTiWiiBii into Z.block(i, m +j) Sum EiTiWijEii into H.block(i, 0) : Sum BiTiWijEii into H.block(m +], 0) } //end i } // end]
11. {solveNormalEquation} Form the lower transpose of the Z matrix and solve the system A = Z-1H. Note that the normal equation system is a symmetric system (e.g., the normal matrix Z is symmetric). Thus, a symmetric system solver can be used. In the case of a symmetric system solver, it may not be necessary to form the lower triangle.

Date Recue/Date Received 2021-02-10 12. {updateParameters} Update all the parameters by extracting the corrections from the A matrix as in Equations 51 and 52.
13. If (p) # 0 compare with the previous root mean square (RMS) of residuals and check for convergence. The convergence condition can be IR(P) _ R(P-1-)1 < 6 Equation 103 14. {computePostConvergenceResiduals + checkForBlunders} If convergence has been reached, perform automatic blunder editing. If convergence has not been reached, increment the iteration index (p) + 1) Equation 104 and go to step 7.
[00201] What follows is pseudocode for the operation 818 for building the reduced normal equations system, computing corrections to ground point positions and performing error propagation via extraction of data from the reduced normal matrix.
This portion of the pseudocode includes ground point coordinate folding.
[00202] As in the full normal solution provided in the previous pseudocode, the same per-image elements are cached in a workspace object and updated with each iteration. The algorithm for the reduced solution can be broken into two major portions:
priming and folding. Priming involves storing of weights and the contributions along the diagonal of the full normal equation matrix (and corresponding data for the right hand column vector H). This corresponds to the 2 portion of Z. Thus, priming involved formation of the minuends of Equation 75 and Equation 76. Folding can include incorporation of the subtrahends of the aforementioned Equations.
[00203] To provide an efficient implementation, a ground point workspace can be created. The workspace can include the following elements: 2j,iii, 2r1-. These things are indexed by ground point for the ground point workspace. The technique can begin by setting the non-linear least squares iteration index (p) to zero.
1. {initializeData} Store nominal initial values for the solved parameters for each image i E tO, === , m ¨ 1). The initial values for Vi( ) and 0[0) can also be set along with eio) as in Equation 98.
2. {initializeData} Compute a. Initial ground point coordinates via Equation 11 for each ground point j E
tO, === , n ¨ 1). These form ).

Date Recue/Date Received 2021-02-10 b. Initial image cache data as described above.
3. {outputIterationReport} Output initial iteration report (input parameters, initial ground point, and ground point observation coordinates).
4. {initializeMatrices} Block partition the reduced normal matrix Minto m subblocks, each of size 6 x 6. Block partition the reduced column vector C
similarly 5. {computeAndStoreDiscrepancies} For each ground point index j E
{0, ===, n ¨ 1) and each observation of that ground point i E 0i, compute the discrepancy vector Eij given in Equation 29 as:
for] E {0, === , n ¨
Fetch the most recent ground point position for i E I and observation index b E
a) Retrieve image-cached values for Ti and 17i b) Retrieve ground point observation 1'j =
c) Apply the observation equation to obtain the projected value for 17, from Equation 99.
d) Compute and store the discrepancy vector for observation b as in Equation 100 end for i Wend for j 6. Compute the standard error of unit weight.as in Equation 57.
7. {buildNormalMatrices}
a. Zero the reduced matrices M and C
b. { initialize Weights} Initialize M (and likewise C) with the block weight matrices on the diagonal. This involves setting the blocks of M to the subblocks of W
and setting the subblocks (subrows) of C to ¨we.
o. {sumInPartialsAndFoldGPs} Form the main diagonal and ground point matrices 2i by iterating over ground points. Perform folding for ground point j for] E {0, ===, n ¨

Date Recue/Date Received 2021-02-10 PRIMING:
Store VV./ into 2./ of GPWS
Store ¨VV./C./ into fif of GPWS
for i E I (where I is set of image indexes upon which GP j is an observation) Form partial derivatives:
Build Pi./ as in Equation 101 Build B11 as in Equation 102 Retrieve discrepancy vector Eif as computed in Equation 100.
Retrieve observation weight matrix Sum in contribution of GP j's obs in image i within M:
Sum EiTiVT7i./Pii into M.block(i,i) Sum EiTiVIVii into C.block(1,0) Sum in i's contribution to 2./ and Sum into 2./
Sum biTilii7Eif into end i Invert 2./ and store into GPWS as 2I1 FOLDING INTO M (note : iteration loop over] is still valid) for r E
Form Zri = EiTitiViiBii as in Equations 69 and 47 for c E I c>r Form 2,.T./
Sum in ¨2ri2I12cTi into M.block(r,c).
Wend c Sum in ¨2r./2111:1./ into C. block(r, 0).
Wend r III end j Date Recue/Date Received 2021-02-10 8. Complete the lower diagonal entries of M and solve A = C. As in the full normal equation solution, note that M is symmetric and thus a symmetric system solver is in order.
9. First use pseudocode provided below to compute corrections to ground points. Then update all the parameters from the A vector.
10. If (p) # 0 compare with the previous RMS of residuals and check for convergence. The convergence condition is IR(P) _ R(p-1)I < E Equation 105 11. {computePostConvergenceResiduals + checkForBlunders} If convergence has been reached, perform automatic blunder editing as detailed elsewhere. After there are no more blunders, proceed with error propagation {propagateErrors}. If convergence has not been reached, increment the iteration index as in Equation 104 and go to step 5.
[00204] After the solution vector A is obtained, unfolding the ground point corrections is a matter of employing Equation 80, replicated here for reference:
A. = - tE.Eli 2T.A.1 Equation 80 1 tj t 3x1 {unfoldGroundPoints}
for] E 0, = = = , n ¨
Retrieve 21' and from GPWS
Store 1.-1 into a new subtrahend matrix S (i.e. initialize S to for r E
Form Zri SUM ¨2rTiAr into S.
} I I end r Compute Ai =
Wend j [00205] The general cross error covariance between ground point indexes r and c can obtained by evaluation of Equation 94.
{relativeGroundPointCov}
Retrieve 2,71- and 2 c-1- from ground point workspace.

Date Recue/Date Received 2021-02-10 Obtain indexing sets /r and /, (image indexes of ground points r and c) Allocate matrix lr,c and initialize to zero. This is the output of this function.
3x3 Allocate matrix P and initialize to zero 3x7 for t E
{
Allocate matrix Q and initialize to zero 3x7 for s E ir {
Form 2,,,, and 2tcfrom Equation 47. (Note that for the reduced system, N =
2) Extract kt=i".getBlock( s, t), where is as defined in Equation 93 Set Q <¨ Q + nti'st Wend s Set P <¨ P + Q2, Wend t Compute 1r,c = 2; 1P2c71 if( r = c) {
Set fr,c, <¨ [0-(P)]22c71 fr,c, I
return rx //end {relativeGroundPointCov}
[00206] The full 3n x 3n ground covariance matrix 1' may be obtained by 3n x 3n invoking the method for r E f1,2, === , n} and for c E fr, r + 1, = = = , nj.
Note that the indexing for c starts with r since the full ground covariance matrix is symmetric (i.e., build the upper triangle of I and "reflect about the diagonal" to obtain the lower 3n )< 3n symmetric portion).
[00207] What follows regards how to perform operation 1014. The operation 1014 proceeds given the outputs of the MLSE techniques discussed. The compensation applies the inverse of the observation equation, accommodating the various relative Date Recue/Date Received 2021-02-10 frame offsets to arrive at compensated world space coordinates from misregistered world space coordinates Vw [00208] = Compensate,(Vreg) Equation 106 reg [00209] The motivation for providing the inputs and outputs in world space coordinates can be that is the native space of the inputs and desired space of the outputs for each element of each image's point cloud.
[00210] For an arbitrary misregistered vector V,reg on image i, the compensation formula can be performed as in Equation 107 ¨ TiT (VmWisreg T7i(0)R V : ) T71( p) W
reg VR
V (1+ Si(P)) Equation 107 [00211] where T, is constructed from the solution vector 9(,P) and the other symbols in Equation 107 are defined elsewhere. Note that the values for ¨17( )R ¨
and 47,(P) + and TT can be precomputed on a per image basis when applying Equation 107 for a time-efficient implementation.
[00212] FIG. 11 illustrates, by way of example, a diagram of an embodiment of a method 1100 for 3D point set generation and registration. The method 1100 as illustrated includes capturing, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of the object, at operation 1102;
registering the overlapping subsections to each other, at operation 1104; and geo-locating the registered overlapping subsections, at operation 1106.
[00213] The method 1100 can include capturing, by a UV of the UVs, a first overhead image of a starting geo-location at which the image data is captured and wherein geo-locating the overlapping subsections includes correlating the first overhead image with a second overhead image for which geo-location is known. The method 1100 can include, wherein the second overhead image is a satellite image. The method 1100 can include, wherein geo-locating the registered overlapping subsection includes determining a normalized cross correlation of image chips of the first overhead image and the second overhead image.
[00214] The method 1100 can include receiving, from an operator of a UV of the UVs, a starting geo-location, and a heading of the UV. The method 1100 can include wherein geo-locating the registered overlapping subsections is performed based on the Date Recue/Date Received 2021-02-10 starting geo-location and the heading. The method 1100 can include performing, by a UV of the UVs, a light detection and ranging (LIDAR) scan to generate LIDAR
scan data. The method 1100 can include wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.
[00215] The method 1100 can include associating, by the UV, geo-location data of the UV with image data generated by the UV. The method 1100 can include, wherein geo-locating the registered overlapping subsections occurs based on the geo-location data. The method 1100 can include generating a first three-dimensional (3D) point set based on the geo-located registered overlapping subsections. The method 1100 can include registering the first 3D point set to a second 3D point set to generate a merged 3D point set. The method 1100 can include, wherein registering the first 3D
point set to the second 3D point set includes scaling, rotating, and translating one or more of the first and second 3D point sets using a least squares estimate bundle adjustment based on tie points between the first and second 3D point sets.
[00216] FIG. 12 illustrates, by way of example, a block diagram of an embodiment of a machine in the example form of a computer system 1200 within which instructions, for causing the machine to perform any one or more of the methodologies discussed herein, may be executed. In alternative embodiments, the machine operates as a standalone device or may be connected (e.g., networked) to other machines. In a networked deployment, the machine may operate in the capacity of a server or a client machine in server-client network environment, or as a peer machine in a peer-to-peer (or distributed) network environment.
The machine may be a personal computer (PC), a tablet PC, a set-top box (STB), a Personal Digital Assistant (PDA), a cellular telephone, a web appliance, a network router, switch or bridge, or any machine capable of executing instructions (sequential or otherwise) that specify actions to be taken by that machine. Further, while only a single machine is illustrated, the term -machine" shall also be taken to include any collection of machines that individually or jointly execute a set (or multiple sets) of instructions to perform any one or more of the methodologies discussed herein.
[00217] The example computer system 1200 includes a processor 1202 (e.g., processing circuitry 118, such as can include a central processing unit (CPU), a graphics processing unit (GPU), field programmable gate array (FPGA), other circuitry, such as one or more transistors, resistors, capacitors, inductors, diodes, regulators, switches, multiplexers, power devices, logic gates (e.g., AND, OR, XOR, negate, etc.), buffers, memory devices, or the like, or a combination thereof), a main memory 1204 and a static memory 1206, which Date Recue/Date Received 2021-02-10 communicate with each other via a bus 1208. The computer system 1200 may further include a display device 1210 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)). The computer system 1200 also includes an alphanumeric input device 1212 (e.g., a keyboard), a user interface (UI) navigation device 1214 (e.g., a mouse), a disk drive unit 1216, a signal generation device 1218 (e.g., a speaker), a network interface device 1220, and radios 1230 such as Bluetooth, WWAN, WLAN, and NFC, permitting the application of security controls on such protocols.
[00218] The disk drive unit 1216 includes a machine-readable medium 1222 on which is stored one or more sets of instructions and data structures (e.g., software) 1224 embodying or utilized by any one or more of the methodologies or functions described herein. The instructions 1224 may also reside, completely or at least partially, within the main memory 1204 and/or within the processor 1202 during execution thereof by the computer system 1200, the main memory 1204 and the processor 1202 also constituting machine-readable media.
[00219] While the machine-readable medium 1222 is shown in an example embodiment to be a single medium, the term "machine-readable medium" may include a single medium or multiple media (e.g., a centralized or distributed database, and/or associated caches and servers) that store the one or more instructions or data structures. The term "machine-readable medium" shall also be taken to include any tangible medium that is capable of storing, encoding or carrying instructions for execution by the machine and that cause the machine to perform any one or more of the methodologies of the present invention, or that is capable of storing, encoding or carrying data structures utilized by or associated with such instructions. The term "machine-readable medium" shall accordingly be taken to include, but not be limited to, solid-state memories, and optical and magnetic media. Specific examples of machine-readable media include non-volatile memory, including by way of example semiconductor memory devices, e.g., Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks.
[00220] The instructions 1224 may further be transmitted or received over a communications network 1226 using a transmission medium. The instructions 1224 may be transmitted using the network interface device 1220 and any one of a number of well-known transfer protocols (e.g., HTTP). Examples of communication networks include a local area network (-LAN"), a wide area network (-WAN"), the Internet, mobile telephone networks, Date Recue/Date Received 2021-02-10 Plain Old Telephone (POTS) networks, and wireless data networks (e.g.. WiFi and WiMax networks). The term "transmission medium" shall be taken to include any intangible medium that is capable of storing, encoding or carrying instructions for execution by the machine, and includes digital or analog communications signals or other intangible media to facilitate communication of such software.
[00221] Additional Notes and Examples [00222] Example 1 includes a method for generating a three-dimensional (3D) point cloud of an object, the method comprising capturing, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of the object, registering the overlapping subsections to each other, and geo-locating the registered overlapping subsections.
[00223] In Example 2, Example 1 can further include capturing, by a UV of the UVs, a first overhead image of a starting geo-location at which the image data is captured and wherein geo-locating the overlapping subsections includes correlating the first overhead image with a second overhead image for which geo-location is known.
[00224] In Example 3, Example 2 can further include, wherein the second overhead image is a satellite image.
[00225] In Example 4, Example 3 can further include, wherein geo-locating the registered overlapping subsection includes determining a normalized cross correlation of image chips of the first overhead image and the second overhead image.
[00226] In Example 5 at least one of Examples 1-4 can further include receiving, from an operator of a UV of the UVs, a starting geo-location, and a heading of the UV, and wherein geo-locating the registered overlapping subsections is performed based on the starting geo-location and the heading.
[00227] In Example 6, at least one of Examples 2-5 can further include performing, by a UV of the UVs, a light detection and ranging (LIDAR) scan to generate LIDAR scan data, and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.
[00228] In Example 7, at least one of Examples 1-6 can further include associating, by the UV, geo-location data of the UV with image data generated by the UV, and wherein geo-locating the registered overlapping subsections occurs based on the geo-location data.
[00229] In Example 8, at least one of Examples 1-7 can further include generating a first three-dimensional (3D) point set based on the geo-located registered overlapping Date Recue/Date Received 2021-02-10 subsections and registering the first 3D point set to a second 3D point set to generate a merged 3D point set.
[00230] In Example 9, Example 8 can further include, wherein registering the first 3D
point set to the second 3D point set includes scaling, rotating, and translating one or more of the first and second 3D point sets using a least squares estimate bundle adjustment based on tie points between the first and second 3D point sets.
[00231] Example 10 includes a system comprising unmanned vehicles configured to capture image data representative of respective overlapping subsections of an object, and processing circuitry configured to register the overlapping subsections to each other, and geo-locate the registered overlapping subsections.
[00232] In Example 11, Example 10 can further include, wherein a UV of the UVs is further configured to capture a first overhead image of a starting geo-location at which the image data is captured and wherein geo-locating the overlapping subsections includes correlating the first overhead image with a second overhead image for which geo-location is known.
[00233] In Example 12, Example 11 can further include, wherein the second overhead image is a satellite image.
[00234] In Example 13, Example 12 can further include, wherein geo-locating the registered overlapping subsection includes determining a normalized cross correlation of image chips of the first overhead image and the second overhead image.
[00235] In Example 14, at least one of Examples 10-13 can further include, wherein the processing circuitry is further configured to receive, from an operator of a UV of the UVs, a starting geo-location and a heading of the UAV, and wherein geo-locating the registered overlapping subsections is performed based on the starting geo-location and the heading.
[00236] In Example 15, at least one of Examples 11-14 can further include, wherein a UV of the UVs is further configured to perform a light detection and ranging (LIDAR) scan to generate LIDAR scan data; and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.
[00237] Example 16 includes a (e.g., non-transitory) machine-readable medium including instructions that, when executed by a machine, cause the machine to perform operations comprising receiving, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of an object, registering the overlapping subsections to each other, and geo-locating the registered overlapping subsections.

Date Recue/Date Received 2021-02-10 [00238] In Example 17, Example 16 can further include, wherein the operations further comprise receiving, by the UV, geo-location data of the UV associated with the image data generated by the UV, and wherein geo-locating the registered overlapping subsections occurs based on the geo-location data.
[00239] In Example 18, at least one of Examples 16-17 can further include, wherein the operations further comprise generating a first three-dimensional (3D) point set based on the geo-located registered overlapping subsections and registering the first 3D point set to a second 3D point set to generate a merged 3D point set.
[00240] In Example 19, Example 18 can further include, wherein registering the first 3D point set to the second 3D point set includes scaling, rotating, and translating one or more of the first and second 3D point sets using a least squares estimate bundle adjustment based on tie points between the first and second 3D point sets.
[00241] In Example 20, at least one of Examples 16-19 can further include, wherein the operations further comprise receiving light detection and ranging (LIDAR) scan data of the object from a UV of the UVs; and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.
[00242] Although an embodiment has been described with reference to specific example embodiments, it will be evident that various modifications and changes may be made to these embodiments without departing from the broader spirit and scope of the invention. Accordingly, the specification and drawings are to be regarded in an illustrative rather than a restrictive sense. The accompanying drawings that form a part hereof, show by way of illustration, and not of limitation, specific embodiments in which the subject matter may be practiced. The embodiments illustrated are described in sufficient detail to enable those skilled in the art to practice the teachings disclosed herein. Other embodiments may be utilized and derived therefrom, such that structural and logical substitutions and changes may be made without departing from the scope of this disclosure. This Detailed Description, therefore, is not to be taken in a limiting sense, and the scope of various embodiments is defined only by the appended claims, along with the full range of equivalents to which such claims are entitled.

Date Recue/Date Received 2021-02-10

Claims (20)

What is claimed is:

1. A method for generating a three-dimensional (3D) point cloud of an object, the method comprising:
capturing, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of the object;
registering the overlapping subsections to each other; and geo-locating the registered overlapping subsections.

2. The method of claim 1, further comprising:
capturing, by a UV of the UVs, a first overhead image of a starting geo-location at which the image data is captured; and wherein geo-locating the overlapping subsections includes correlating the first overhead image with a second overhead image for which geo-location is known.

3. The method of claim 2, wherein the second overhead image is a satellite image generated by a satellite.

4. The method of claim 3, wherein geo-locating the registered overlapping subsection includes determining a normalized cross correlation between image chips of the first overhead image and image chips of the second overhead image.

5. The method of claim 1, further comprising:
receiving, from an operator of a UV of the UVs, a starting geo-location, and a heading of the UV; and wherein geo-locating the registered overlapping subsections is performed based on the starting geo-location and the heading.

Date Recue/Date Received 2021-02-10

6. The method of claim 2, further comprising:
performing, by a UV of the UVs, a light detection and ranging (LIDAR) scan to generate LIDAR scan data; and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.

7. The method of claim 1, further comprising:
associating, by the UV, geo-location data of the UV with image data generated by the UV; and wherein geo-locating the registered overlapping subsections occurs based on the geo-location data.

8. The method of claim 1, further comprising:
generating a first three-dimensional (3D) point set based on the geo-located registered overlapping subsections; and registering the first 3D point set to a second 3D point set to generate a merged 3D
point set.

9. The method of claim 8, wherein registering the first 3D point set to the second 3D
point set includes scaling, rotating, and translating one or more of the first and second 3D
point sets using a least squares estimate bundle adjustment based on tie points between the first and second 3D point sets.

10. A system comprising:
unmanned vehicles configured to capture image data representative of respective overlapping subsections of an object; and processing circuitry configured to:
register the overlapping subsections to each other; and geo-locate the registered overlapping subsections.

11. The system of claim 10, wherein a UV of the UVs is further configured to capture a first overhead image of a starting geo-location at which the image data is captured and wherein geo-locating the overlapping subsections includes correlating the first overhead image with a second overhead image for which geo-location is known.
Date Recue/Date Received 2021-02-10

12. The system of claim 11, wherein the second overhead image is a satellite image.

13. The system of claim 12, wherein geo-locating the registered overlapping subsection includes determining a normalized cross correlation of image chips of the first overhead image and the second overhead image.

14. The system of claim 10, wherein the processing circuitry is further configured to:
receive, from an operator of a UV of the UVs, a starting geo-location, and a heading of the UAV; and wherein geo-locating the registered overlapping subsections is performed based on the starting geo-location and the heading.

15. The system of claim 11, wherein a UV of the UVs is further configured to perform a light detection and ranging (LIDAR) scan to generate LIDAR scan data; and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.

16. A machine-readable medium including instructions that, when executed by a machine, cause the machine to perform operations comprising:
receiving, by unmanned vehicles (UVs), image data representative of respective overlapping subsections of an object;
registering the overlapping subsections to each other; and geo-locating the registered overlapping subsections.

17. The machine-readable medium of claim 16, wherein the operations further comprise:
receiving, by the UV, geo-location data of the UV associated with the image data generated by the UV; and wherein geo-locating the registered overlapping subsections occurs based on the geo-location data.

18. The machine-readable medium of claim 16, wherein the operations further comprise:
generating a first three-dimensional (3D) point set based on the geo-located registered overlapping subsections; and Date Recue/Date Received 2021-02-10 registering the first 3D point set to a second 3D point set to generate a merged 3D
point set.

19. The machine-readable medium of claim 18, wherein registering the first 3D point set to the second 3D point set includes scaling, rotating, and translating one or more of the first and second 3D point sets using a least squares estimate bundle adjustment based on tie points between the first and second 3D point sets.

20. The machine-readable medium of claim 16, wherein the operations further comprise receiving light detection and ranging (LIDAR) scan data of the object from a UV of the UVs;
and wherein geo-locating the registered overlapping subsections includes correlating the first overhead image with the LIDAR scan data.

Date Recue/Date Received 2021-02-10