Orthophoto image mosaic method using adaptive cost A* algorithm

Abstract

This invention discloses an orthophoto mosaicking method using an adaptive cost A* algorithm. The method includes obtaining initial mosaicking lines between pairs of orthophotos to be mosaicked to construct an initial mosaicking line network; optimizing the initial mosaicking line network using an improved A* algorithm, i.e., using an adaptive threshold algorithm to determine traversable areas, and constructing a novel cost function to obtain the improved A* algorithm; then using the improved A* algorithm to optimize each initial mosaicking line in the initial mosaicking line network to obtain optimized mosaicking lines, thereby completing the mosaicking between orthophotos; the adaptive threshold algorithm uses the average grayscale value of the overlapping area of ​​adjacent images as the global threshold, and the smaller of the average grayscale values ​​of N*N blocks in the left and right images where the current search point is located as the local threshold. Based on the maximum grayscale difference cost of the current search point, it adaptively selects either the global threshold or the local threshold as the adaptive threshold to determine whether the current search point is a traversable area.

Description

Technical Field

[0001] This invention belongs to the technical field of path planning, specifically relating to an orthophoto mosaicking method using the adaptive cost A* algorithm. Background Technology

[0002] Orthophotos (DOMs) provide high-precision geographic information and are widely used in cartography, urban planning, environmental monitoring, and other fields. However, the image range of a single DOM is limited, and large-scale DOMs are usually mosaicked. Due to the limitations of DOM generation accuracy, directly superimposed and stitched orthophotos will show obvious misalignment of the same features at the edges. Inconsistencies in geometric and radiometric phase between adjacent images will lead to the mosaicking line problem, which affects the quality of large-scale orthophotos. To solve this problem, research mainly focuses on optimizing the mosaicking lines between adjacent images to minimize radiometric and geometric differences. The basic steps can be divided into two stages: the first stage extracts grayscale, texture, and target information from the overlapping areas of the orthophotos and abstracts the images into a cost map based on these features; the second stage optimizes the energy function based on the cost map, searches for mosaicking lines, and generates a seamless large-scale DOM. From a radiometric perspective, the goal is to reduce color inconsistencies; geometrically, the aim is to eliminate the misalignment of features near the mosaicking lines.

[0003] However, many existing methods for large-scale image mosaicking often perform well only in specific regions when dealing with different land cover types. Specific algorithms are only applicable to densely built-up urban areas or sparsely textured mountainous and forested areas, with limited adaptability to complex terrain and diverse land cover. In addition, guided methods based on external data (such as road and building vectors) rely on the availability and accuracy of auxiliary data, limiting their application in areas with scarce or low-quality data. These factors result in a lack of general and robust algorithms for large-scale image mosaicking tasks, thereby increasing the cost of manual intervention and operation.

[0004] Therefore, there is an urgent need for a robust method for large-scale remote sensing image mosaicking that can effectively handle complex scenes containing various land features such as mountains, rivers, and houses. Summary of the Invention

[0005] This invention proposes an orthophoto mosaicking method based on the adaptive cost A* algorithm to address the challenge of mosaicking line detection in large-scale remote sensing image mosaicking. This method aims to significantly improve the robustness and computational efficiency of image mosaicking. Through adaptive threshold judgment and heuristic path search mechanisms, it effectively optimizes mosaicking lines, reducing the impact of geometric distortion and radiometric inconsistencies. This makes the algorithm applicable to mosaicking line extraction in remote sensing images containing both densely built-up urban areas and sparsely built-up mountainous regions.

[0006] This invention can be achieved through the following technical solutions:

[0007] An orthophoto mosaicking method using the adaptive cost A* algorithm includes the following steps:

[0008] Step 1: Obtain the initial tessellation lines between each pair of orthophoto images to be tessellated, in order to construct the initial tessellation line network;

[0009] Step 2: Optimize the initial tessellation network using an improved A* algorithm.

[0010] An adaptive threshold algorithm is used to determine the passable area and a new cost function is constructed to obtain an improved A* algorithm. Then, the improved A* algorithm is used to optimize the initial mosaic lines in the initial mosaic line network one by one to obtain the optimized mosaic lines, thereby completing the mosaicking between orthophoto images.

[0011] The global threshold is the average gray value of the overlapping area of ​​adjacent images, and the local threshold is the smaller of the average gray values ​​of the N*N blocks in the left and right images where the current search point is located. The adaptive threshold algorithm adaptively selects either the global threshold or the local threshold as the adaptive threshold and compares it with the maximum gray value difference cost of the current search point to determine whether the current search point is a passable area.

[0012] The novel cost function calculates the grayscale cost based on the maximum grayscale difference cost of the current search point, and improves the cost function by incorporating distance control costs.

[0013] Furthermore, when using the adaptive threshold algorithm to determine passable areas, the adaptive threshold T is first calculated using the following formula:

[0014]

[0015]

[0016] in, The global threshold is the average grayscale value of the overlapping region between adjacent images. Representing the current image point The average grayscale value of the N*N blocks in the adjacent left and right images is calculated using formula #2. Indicates the current search point grayscale value;

[0017] choose The smaller one is used as the local threshold. If the calculated local threshold is less than the global threshold, the local threshold is used as the adaptive threshold T; otherwise, the global threshold is used as the adaptive threshold T.

[0018] Then, using the following formula, calculate the maximum grayscale difference cost corresponding to the current search point.

[0019]

[0020]

[0021] in, These represent the current search points in adjacent left and right images, respectively. grayscale value, These represent the current search points in adjacent left and right images, respectively. The maximum gray-level gradient in the eight directional neighborhoods: top left, top, top right, right, bottom right, bottom, bottom left, and left.

[0022] Finally, if the current search point Corresponding maximum grayscale difference cost If the value is greater than the adaptive threshold T, the pixel is determined to be a non-passable area; otherwise, the pixel is determined to be a passable area.

[0023] Furthermore, when using the improved A* algorithm for searching, an adaptive threshold algorithm is first used to determine whether the current search point is a passable area. If not, the search proceeds to the next pixel.

[0024] If so, then search the pixels in the eight-directional neighborhood around the current search point. If any pixel satisfying formulas #5 and #6 is found, then obtain the grayscale cost. ,

[0025]

[0026]

[0027] Among them, adaptive threshold and maximum grayscale difference cost All calculations utilize the method described above to calculate any pixel within the eight-directional neighborhood of the current search point. The corresponding value;

[0028] Then, using the following formula, combined with the above grayscale cost... Calculate the new cost function F and perform subsequent search steps.

[0029]

[0030]

[0031] Where D represents distance control cost, G represents the distance from the starting point to the current search point, H represents the estimated distance from the current search point to the destination, w represents the estimated path cost weight, Grad represents the grayscale cost, c represents the grayscale conversion coefficient used to balance the weight settings, and SSIM represents the structural similarity index. Let α represent the Euclidean distance from the current search point p to the starting point, and let α represent a positive integer greater than 1.

[0032] If no pixel satisfying formulas #5 and #6 is found, the search continues for the next pixel until the search is complete and the optimized mosaic line is obtained.

[0033] Furthermore, the distance control cost D is calculated using the following formula.

[0034]

[0035] in, Indicates the coordinates of the current search point. Represents any line segment between two adjacent points on the initial tessellation line. For each line segment... Calculate the current search point respectively The perpendicular distance to the line segment is .

[0036] Furthermore, when using the improved A* algorithm for the search, the start and end points of the initial tessellation line are adjusted using the following method:

[0037] First, for the current search point The adaptive threshold for the current search point is calculated using the calculation method described above. Intermediate grayscale difference cost ,

[0038] Then, perform a condition check: 1. If Then determine the current search point. If it is a non-obstacle point, then it is an obstacle point; otherwise, it is an obstacle point.

[0039] 2. Traverse from the current search point Let be the pixels within the circle's center and radius dN. If all these pixels are obstacle points, then the current search point... Passable depth At least dN, that is ;

[0040] Finally, if both of the above conditions are met simultaneously and This indicates the current search point. These are non-obstacle points with a wide passable area, so no adjustment is needed;

[0041] Otherwise, perform a depth-first search of the surrounding neighborhood until a pixel that satisfies the above two conditions is found, and use that pixel to replace the current search point. Adjustments completed.

[0042] Furthermore, when constructing the initial mosaicking network, the orthophoto image to be mosaicked first needs to be preprocessed, including geometric correction and downsampling.

[0043] The beneficial technical effects of this invention are as follows:

[0044] 1. Adaptive threshold mechanism:

[0045] This invention dynamically calculates an adaptive threshold based on local and global grayscale information of overlapping image regions to accurately distinguish between passable and impassable areas. Compared to fixed threshold methods, this adaptive approach more accurately identifies feature boundaries and texture characteristics, preventing mosaic lines from crossing obvious feature areas such as buildings, roads, and ridges, effectively reducing the possibility of geometric distortion.

[0046] 2. Improved A* heuristic search strategy:

[0047] This method employs an improved A-algorithm, based on an adaptive threshold, and integrates a heuristic function combining structural similarity (SSIM) and distance-controlled cost. This ensures that the search path not only avoids areas with significant geometric and radial differences but also effectively tracks roads, water areas, and open regions with continuous edges or low grayscale differences. Compared to the traditional A-algorithm, the design of this heuristic function significantly reduces the search space and the number of computational nodes, thereby substantially improving computational efficiency.

[0048] 3. Dynamically constructing tessellation networks:

[0049] The algorithm first uses the polygon skeleton (straight skeleton) method to generate an initial tessellation network, providing an initial search framework. Then, it combines the improved A* algorithm with a heuristic search strategy for fine optimization. This method ensures the reasonable layout of the initial path in the overlapping area, avoids the waste of computational resources caused by the unrestricted search range in traditional methods, and improves the robustness of tessellation search.

[0050] 4. Does not rely on additional auxiliary data:

[0051] This invention utilizes only the grayscale information and structural similarity of the image itself, without requiring a digital surface model (DSM), vector data, or other auxiliary data. This makes the algorithm applicable to areas with large-scale and multi-type land cover, especially in urban areas with high-density buildings and mountainous areas with complex terrain. It effectively makes up for the shortcomings of traditional methods in terms of mosaic line accuracy due to the lack of high-precision auxiliary data.

[0052] 5. Image downsampling improves efficiency:

[0053] By performing Gaussian pyramid downsampling on overlapping image regions, the amount of data in the search area is reduced, effectively improving the algorithm's computational efficiency. The downsampling process simultaneously preserves low-frequency information (outline structure) of ground features, reducing noise and detail interference, ensuring the algorithm can quickly and robustly identify key ground features.

[0054] Through these comprehensive designs, the adaptive cost A* algorithm proposed in this invention can maintain a highly stable mosaic line extraction effect in complex terrain areas with significant geometric distortion and radiometric inconsistency. It solves the problem that traditional methods cannot simultaneously balance the accuracy and efficiency of mosaic line extraction in densely built urban areas and sparsely built mountainous areas. It can solve the problem of mosaic line detection in large-scale remote sensing image mosaicking, and significantly improve robustness and computational efficiency. Attached Figure Description

[0055] Figure 1 This is a structural diagram of the method proposed in this invention;

[0056] Figure 2 This is a schematic diagram of the mosaic line construction process of the input image after downsampling according to the present invention;

[0057] Figure 3 This is a flowchart illustrating the generation of the initial mosaic lines for this invention;

[0058] Figure 4 This is a schematic diagram of the A* algorithm search execution process of the present invention;

[0059] Figure 5 This is a schematic diagram illustrating how the impassable area changes accordingly when N is different in this invention;

[0060] Figure 6 This is a comparative schematic diagram showing the endpoints before and after adjustment when optimizing the initial inlay line according to the present invention;

[0061] Figure 7 This is a schematic diagram illustrating the process of adjusting the endpoints during the optimization of the initial inlay line according to the present invention;

[0062] Figure 8 This is a schematic diagram illustrating the effect of structural similarity in this invention;

[0063] Figure 9 This diagram illustrates a comparison of the algorithm execution efficiency when the integrated image is downsampled at the same rate as the cost constraint and when SSIM is calculated additionally.

[0064] Figure 10 This is a schematic diagram of an orthophoto of a certain region and an initial mosaic network in a specific embodiment of the present invention;

[0065] Figure 11This is a schematic diagram of the optimized mosaic line portion in a certain region in a specific embodiment of the present invention;

[0066] Figure 12 This is a schematic diagram showing the results of urban area mosaic line extraction in a specific embodiment of the present invention;

[0067] Figure 13 This is a comparative diagram of different algorithms used in densely built-up urban areas in a specific embodiment of the present invention;

[0068] Figure 14 This is a comparative diagram showing the application of different algorithms to sparsely built lakes and mountainous areas in a specific embodiment of the present invention.

[0069] Figure 15 This is a schematic diagram demonstrating the effect of geometric consistency in a specific embodiment of the present invention. Detailed Implementation

[0070] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings and preferred embodiments.

[0071] like Figure 1 As shown, this invention provides an orthophoto mosaicking method using an adaptive cost A* algorithm. First, initial mosaicking lines between each pair of orthophotos to be mosaicked are obtained to construct an initial mosaicking line network. Then, an improved A* algorithm is used to optimize the initial mosaicking line network. Specifically, an adaptive threshold algorithm is used to determine passable areas, and a novel cost function is constructed to obtain the improved A* algorithm. Finally, the improved A* algorithm is used to optimize any initial mosaicking line in the initial mosaicking line network to obtain optimized mosaicking lines, thereby completing the mosaicking between orthophotos. The orthophoto mosaicking method proposed in this invention avoids reliance on external auxiliary data. Considering that remote sensing images typically have large coverage and high resolution, it must balance computational efficiency and processing performance. This will help achieve high-quality, low-cost image mosaicking under various geographical environments and data conditions, meeting diverse needs in practical applications.

[0072] Specifically as follows:

[0073] Step 1: Obtain the initial tessellation network

[0074] 1. Preprocess the acquired images

[0075] Because high-resolution satellite remote sensing imagery has a large image width, a single image typically contains a large amount of data. Since the A* algorithm's cost calculation is based on the number of pixels searched, the time required for mosaic path planning increases with the amount of image data. To improve the efficiency of mosaic optimization, we preprocess the overlapping areas of the images, including geometric correction and downsampling such as pyramid downsampling, to reduce geometric distortion in large-scale remote sensing images. This reduces the data scale of the path search, laying the foundation for subsequent mosaic extraction and optimization.

[0076] Typically, the downsampling ratio is 1:3 for aerial imagery, and can be increased to 1:5 for wide-swath satellite imagery. During downsampling, we utilize Gaussian convolution kernels to effectively suppress high-frequency information in the remote sensing imagery. High-frequency information usually includes details and noise, while low-frequency information represents the overall structure and contours of the image. By using Gaussian kernels to process the image, the stability of true edge recognition can be enhanced, thereby improving the accuracy and reliability of feature extraction. We downsample the overlapping regions and convert the coordinates of the initial mosaic lines to the downsampled coordinates. After optimizing the mosaic lines on the downsampled overlapping regions, the coordinates of the optimized mosaic line path points are converted back to the coordinates of the original image, and the optimized mosaic lines replace the initial mosaic lines, completing the optimization of the mosaic line network (see Appendix). Figure 2 .

[0077] 2. Generate the initial tessellation network

[0078] By extracting the effective region contours of the orthophoto and further extracting the contours of the overlapping regions between the images, both of which are simple convex polygons, the skeleton bisectors of the polygons in the extracted overlapping regions are calculated as initial mosaic lines. Then, by analyzing the geometric relationships and boundaries of the overlapping regions, a preliminary mosaic line network is generated to provide input for subsequent optimization steps.

[0079] The initial mosaicking network of this invention is constructed using the Straight Skeleton method in the Computational Geometry Algorithm Library (CGAL). The construction process of the polygon skeleton is based on the Wavefront Propagation method. First, a simple polygon is input, and its boundary is regarded as the initial wavefront, which is prepared to move inward uniformly in a manner parallel to the original edges. During the inward propagation of the wavefront, the algorithm continuously detects and processes events such as vertex collisions and edge vanishing, dynamically updating the wavefront state to ensure the correctness of the inward movement. During this process, key geometric events are recorded and connected to generate a complete straight skeleton structure. When the wavefront completely shrinks to the centroid of the polygon, or when all boundaries shrink inward and vanish, the algorithm terminates, and the straight skeleton structure of the polygon is finally obtained. Finally, when multiple contours shrink to a point or a line segment, the skeleton point can be generated. In practical applications, for multiple rectangular bounding boxes with a common overlapping area, we first extract skeleton lines for each pair of bounding boxes. Utilizing the characteristic that these skeleton lines converge at a point within the overlapping area, we obtain an initial mosaicking network through mutual trimming. This provides a robust initial topology for the subsequent stitching process. (See Appendix) Figure 3 .

[0080] Step 2: Optimize the initial tessellation network using an improved A* algorithm.

[0081] When using the improved A* algorithm for searching, the adaptive threshold algorithm is first used to determine whether the current search point is a passable area. If not, the search proceeds to the next pixel.

[0082] If so, then search the pixels in the eight-directional neighborhood around the current search point. If any pixel that satisfies formulas #5 and #6 is found, then obtain the grayscale cost. Then, calculate the new cost function F and continue with the subsequent search steps;

[0083] If no pixel satisfying formulas #5 and #6 is found, the search continues for the next pixel until the search is complete and the optimized mosaic line is obtained.

[0084] It should be further explained that, due to the heuristic nature of the A* algorithm, it does not actually traverse the entire cost graph. Therefore, the determination of impassable regions is not based on a pixel-by-pixel traversal, but rather on dynamically calculating whether the current node is an impassable region during the search process. After setting the start and end points, the algorithm actually calculates the cost while advancing the path search process, as shown in the appendix. Figure 4 .

[0085] (1) Calculate impassable areas

[0086] The adaptive thresholding algorithm uses the average grayscale value of the overlapping region of adjacent images as the global threshold and the smaller of the average grayscale values ​​of the N*N blocks in the left and right images of the current search point as the local threshold. It adaptively selects either the global or local threshold and compares it with the maximum grayscale difference cost of the current search point to determine whether the current search point is a passable area. In this way, the adaptive thresholding algorithm identifies typical features, i.e., obstacles, in the image. Based on the image grayscale values ​​of the overlapping region of adjacent images, texture, gradient, and edge information are extracted to construct the adaptive threshold, providing a high-quality initial cost map representing impassable areas for the subsequent A* algorithm.

[0087] Specifically as follows:

[0088] We describe the mosaicking process as follows: Based on the grayscale information of the overlapping areas of remote sensing images, an initial cost map is constructed. The A* algorithm is then used to create a path within this initial cost map, specifying impassable areas. Since the images used for mosaicking vary in time and angle, to ensure seamless image stitching and avoid feature misalignment, the path must avoid prominent features like houses and follow roads as much as possible. Houses, ridges, and other features can be marked as obstacles and designated as impassable areas. The A* algorithm will avoid these impassable areas during the path search process.

[0089] Using only local thresholds, which calculates the grayscale information of some pixels around the current search point during the search process, can lead to problems where local thresholds cannot effectively identify obstacles in areas with indistinct grayscale changes. Conversely, using only global thresholds, which consider the grayscale of all pixels in the entire overlapping area, can misclassify inconspicuous buildings with small grayscale values ​​as passable areas. Therefore, in the mosaic line optimization process, by setting an appropriate obstacle threshold, we can determine whether the current search point is impassable or passable. If the maximum grayscale difference cost is higher than the obstacle threshold, the pixel is considered impassable; if it is lower, it is considered passable. We employ an adaptive thresholding algorithm, which uses a threshold segmentation strategy to determine the classification of impassable areas during the search process.

[0090] Based on the above considerations, the obstacle threshold, i.e., the adaptive threshold T, is set by combining the local threshold and the global threshold. That is, the adaptive threshold T of the current search pixel is calculated based on the overall gray value of the overlapping area of ​​the adjacent images and the gray value of a small area of ​​N×N blocks around the current point, as shown in Formulas 1 and 2.

[0091]

[0092]

[0093] The average gray level of the overlapping region of adjacent images is: That is, the global threshold, which is the average grayscale value of N*N blocks in the left and right images of the current search point. These are the local thresholds, which are calculated using Formula 2. This represents the grayscale value of the current search point. If the smaller of the local thresholds is less than the global threshold, that local threshold is used as the adaptive threshold; otherwise, the global threshold is used as the adaptive threshold.

[0094] After determining the adaptive threshold T, we need to calculate the grayscale cost of each pixel in the overlapping region of adjacent images. This allows us to compare it with the adaptive threshold T to complete obstacle detection and also facilitates cost estimation using the A* algorithm during subsequent mosaicking path search. To prevent paths from passing through areas with significant intensity, geometric inconsistencies, or buildings within the overlapping region of adjacent images, we should maximize the grayscale cost of these areas, i.e., use the maximum grayscale difference cost to describe the grayscale cost. Grayscale differences between left and right images can detect terrain changes in adjacent images, while grayscale gradations around image pixels can detect areas such as bridges and buildings. Therefore, we can construct the maximum grayscale difference cost for pixels in the overlapping region using grayscale difference and grayscale gradation values. As shown in formulas (3) and (4)

[0095]

[0096]

[0097] in, This represents the grayscale value corresponding to the current search point in the adjacent left and right images. These correspond to the maximum grayscale gradients in the eight directional neighborhoods of the current search point in the adjacent left and right images: top left, top, top right, right, bottom right, bottom, bottom left, and left. If the maximum grayscale difference cost is directly used as the grayscale cost... The calculation function cannot detect continuous building areas with small grayscale differences between the left and right images.

[0098] Furthermore, the adaptive thresholding method based on formulas (1)-(4) can be further improved by using dual thresholds for edge detection and connection to better handle strong and weak edges. Therefore, the grayscale cost in the novel cost function used for subsequent improvements to the A* algorithm... The calculation is as follows:

[0099]

[0100]

[0101] in, This represents the coordinates of any pixel in the eight neighborhoods of the current search point. The maximum grayscale difference cost of any pixel in the neighborhood of the current search point is represented by formulas (3) and (4). and Representing pixels The corresponding gray values ​​in the left and right images It is a pixel. The adaptive threshold is calculated from formulas (1) and (2).

[0102] Cost of maximum grayscale difference All search points with values ​​less than the adaptive threshold T will undergo further auxiliary thresholding, i.e., their original maximum grayscale difference cost will be reassigned. for Gray cost used for subsequent cost function calculation also by Decision. Based on the grayscale cost and grayscale value of any point in the neighborhood, determine whether it is a strong edge. If the maximum grayscale difference cost between any two adjacent points is... Greater than or equal to the adaptive threshold If the pixel has strong edges and no path is found, a higher cost is allocated for subsequent cost function calculations by multiplying the maximum grayscale difference cost corresponding to that pixel by a large constant value of 10000. Similarly, if the grayscale values ​​of any adjacent points in the left and right images are... or Greater than or equal to the adaptive threshold If the grayscale value is twice the value of the edge, the pixel is considered a relatively bright area and a potential secondary edge point with high grayscale intensity, but the edge may be weak. In the absence of a traversable path, a higher cost is allocated for subsequent cost function calculations by multiplying the maximum grayscale difference cost corresponding to this pixel by a large constant value of 2. If neither of these conditions is met, the search point is likely a weak edge, and the original maximum grayscale difference cost is used as the traversal grayscale cost for subsequent cost function calculations. Impassable areas include... Figure 5 As shown, the impassable areas change accordingly when the obtained N is different.

[0103] The improved A* algorithm avoids impassable regions during the search process, i.e. Value ∞ The node, ∞, represents the maximum value allowed in the computer. It should be further noted that... Figure 5The impassable areas shown are obtained by traversing all pixels of the overlapping image. This is only to show the effect of impassable areas. In fact, before the algorithm is executed, it does not completely traverse all pixels of the overlapping image to obtain all impassable areas. Instead, it dynamically determines whether the search point is an impassable area as the A* algorithm searches.

[0104] (2) Endpoint adjustment

[0105] Based on the calculation of impassable regions, the start and end points of the initial mosaicking lines to be optimized can be adjusted. This is because the intersection points of the mosaicking lines in the initial mosaicking network are not always perfect, such as... Figure 6 (a) The intersection of the three initial mosaic lines marked by the blue dashed line box, since the endpoint and start point may be located in obvious building areas, such as Figure 6 As shown in (b), if vertices are located on buildings, the mosaic lines will inevitably pass through them, resulting in geometric misalignment in the image. Therefore, it is crucial to move the points from the rooftop to the ground near the buildings. These points should be confined to a defined area, allowing them to move freely without altering the topology of the initial mosaic line network. To this end, we propose an endpoint adjustment method that uses the obstacle state of the endpoints and their traversable depth as criteria. Adjustment is limited to the overlapping areas of adjacent images. A schematic diagram of the endpoint adjustment is shown below. Figure 6 (c).

[0106] The specific steps are as follows: For the current search point First, the corresponding adaptive thresholds are calculated using formulas (1) and (2). Calculate the grayscale cost using formulas (3)(4)(5)(6) ,

[0107] Then two conditions need to be checked: 1. Check the current search point. Is it an obstacle point? If, then the point is a non-obstacle point, if If , then it is an obstacle point;

[0108] 2. Traverse from the current search point Let be the number of pixels within a circle with center dN and radius dN. If there are no impassable points around a given point, then the passable depth of that point is... At least dN, that is ;

[0109] Therefore, the passable depth Its function is to determine whether there is a path from the current search point to the destination. If not, it means that the current search point is in an area closed by obstacles, such as the roof area, and the current point needs to be adjusted to go outside the impassable area.

[0110] Finally, if both conditions are met and This indicates the current search point. If the point is a non-obstacle point and has a wide passable area, then the current search point is... No adjustment is needed; otherwise, perform a depth-first traversal of the surrounding neighborhood to find new pixels that satisfy the above two conditions, and set the current search point as the new pixel. Adjust to the new pixel, where, r and c represents the row and column coordinate offsets between the new pixel and the current search point. The specific process is as follows: Figure 7 As shown.

[0111] (2) Calculation of a new cost function

[0112] Based on the initial cost map, an improved A* algorithm is used to optimize the initial mosaicking lines. By designing a reasonable heuristic function, the number of search nodes is significantly reduced, improving computational efficiency. The optimized mosaicking lines can effectively avoid high-gradient areas such as building edges and abrupt changes in terrain features, while satisfying geometric and radiometric consistency, ensuring high-quality image mosaicking. Finally, the optimized mosaicking line network is applied to the image mosaicking process for global mosaicking. The generated mosaicking lines are used to create an image mask, ultimately generating a seamless, high-quality, large-scale orthophoto.

[0113] As shown in Equation 8, the cost function of the standard A* algorithm defines F as the total cost from the starting point to the end point through the current position, G as the distance from the starting point to the current point, and H as the estimated distance from the current point to the end point, as shown in Equation 8.

[0114]

[0115] When the A* search algorithm uses a heuristic cost function for shortest path search, the estimated path cost distance H of the current point should be as close as possible to the actual distance from the current point to the destination to make the A* algorithm more accurate. If the H value is based on the Manhattan distance, finding a path with a length less than or equal to that distance is considered a success. However, if the diagonal distance is used, the goal is to find a path with a length greater than or equal to the diagonal distance that is as short as possible. In this method, the search direction includes the eight neighbors of the current point, including both straight and diagonal movements. Euclidean distance is used as the estimated path cost. To improve computational efficiency, integer units are used, with a cost of 10 for a straight movement and 14 for a diagonal movement.

[0116] For orthophoto mosaicking, the mosaicking lines should represent the shortest paths with minimal geometric and radiometric differences. After radiometric and color correction, only geometric inconsistencies need to be considered. To reduce geometric differences, the mosaicking lines should bypass regions with significant grayscale differences and gradations between adjacent images. By controlling the grayscale cost of the search points to avoid such regions, the mosaicking lines can be further optimized by adjusting the initial mosaicking line network using multiple point-to-point distance constraints. Therefore, we can modify the cost function of the A* algorithm and propose a novel cost function, as shown in the following formula:

[0117]

[0118] Where w represents the estimated path cost weight; Grad represents the grayscale cost. When the searched node is determined to be an impassable area, the value of Grad is calculated as ∞, that is, the maximum value allowed in the computer is used in the calculation; c represents the grayscale conversion coefficient used to balance the weight setting, to ensure that the calculation of grayscale cost change is numerically close to the distance unit. In this method, c=30; SSIM represents the structural similarity index; D represents the distance control cost of the initial mosaic line to be optimized, see formula (17).

[0119] During pathfinding, the A* algorithm performs better when the estimated distance H is closer to the actual distance from the current search point to the endpoint. Since obstacles need to be avoided during pathfinding, directly using the Euclidean distance from the current point to the endpoint as the estimated cost distance deviates from the actual search path. The cost of the search path (i.e., the actual path from the current search point to the starting point) can be used as a reference for cost estimation. This method uses the ratio between the actual search path cost from the current point to the starting point and the Euclidean distance from the current point to the starting point as a reference weight for the estimated distance. As the actual search distance from the current point to the starting point increases, the estimated distance cost also increases. The weight is increased by applying a power operation to the reference weight. The formula for calculating the estimated path cost weight w is as follows:

[0120]

[0121] in, This represents the Euclidean distance from the current search point P to the starting point. α is a positive integer greater than 1, intended to increase the derivative of the weight value. It should not be too large; in this method, it is set to 3.

[0122] In the adaptive gray-scale thresholding method for detecting impassable areas, after endpoint adjustment, the mosaic lines can avoid major building outlines and ridgelines. However, this method cannot detect all feature differences. In areas where feature outline information is unclear or gray-scale changes are weak, the adaptive thresholding only detects approximate feature outlines. (See...) Figure 8In (a) and (b), we aim for the mosaic lines to avoid areas with minimal changes in land features within the overlapping regions. To this end, we introduce the Structural Similarity Index (SSIM) as an additional measure of structural similarity to more effectively detect subtle differences. After calculating the SSIM, the mosaic lines will avoid areas with low structural similarity, thus avoiding even smaller differences within the overlapping regions, as shown in [see...]. Figure 8 (c) During the A* algorithm search, a small area of ​​M×M window size is extracted from the two images in the overlapping region of this pixel, and the SSIM of this node is calculated. If the SSIM calculation result is large, it means that the extracted small area images are more similar, and the mosaic line is more likely to pass through this region; if the SSIM calculation result is small, it means that the extracted small area images are not similar, and the mosaic line needs to avoid this region.

[0123] The Structural Similarity Index (SSIM) can be calculated using the "Structural Similarity (SSIM) Index Algorithm" proposed by Wang et al. in 2004. The calculation process includes calculating the mean, standard deviation and covariance between images within a local window, and combining these statistics to obtain a comprehensive similarity index.

[0124] By comparing the brightness of pixel values ​​within the window Contrast and structural similarity SSIM measures the local similarity of images within an M×M window, where x and y represent neighboring images. Typically, M is set to 3. To improve algorithm efficiency, we usually downsample the image; however, SSIM is calculated on the original resolution image. If the image has been downsampled at a 1:3 ratio, the calculation will be based on pixels within a 9×9 area centered at the search point. The formula for SSIM is as follows:

[0125]

[0126] in, , , , and It is the average brightness value within an M×M window. and It is the standard deviation within an M×M window, while It is the covariance within an M×M window; , and It is a constant used to avoid instability when the denominator is very small, and is usually 1. , and Where L is the dynamic range of pixel values, typically [0, 255]; α, β, and γ are parameters that adjust the relative importance of the three components, usually set to α=β=γ=1. The calculation formula is as follows:

[0127]

[0128]

[0129]

[0130]

[0131]

[0132] in, This represents the downsampling factor.

[0133] SSIM can take values ​​in the range [0, 1], where 1 indicates that the two images being compared are the same. Therefore, 2-SSIM takes values ​​in the range [1, 2], which means that the smaller the value, the more similar the two images are and the lower the cost.

[0134] In terms of efficiency, due to the heuristic search nature of the A* algorithm, the cost calculation is also dynamically calculated as the A* algorithm's path search progresses, thus theoretically possessing high efficiency. Furthermore, we evaluated the algorithm's execution efficiency when the composite image's downsampling magnification was increased and when SSIM was additionally calculated as a cost constraint; see [link to relevant documentation]. Figure 9 It can be observed that the main execution efficiency of the algorithm depends on the number of pixels in the overlapping area, and the execution time of the algorithm does not increase significantly after the additional introduction of SSIM.

[0135] During pathfinding using the A* algorithm, a distance control cost is introduced to manage the distance to the initial tessellation line, ensuring that pathpoints are as close as possible to the initial tessellation line and thus preventing excessive deviation of the tessellation line from the initial position. The distance control cost D is defined as the minimum straight-line distance from the current search point to the initial tessellation line, calculated using Euclidean distance.

[0136]

[0137] in, Indicates the coordinates of the current search point. This represents any line segment between adjacent points on an initial tessellation line. An initial tessellation line can be composed of several... Composition. For each line segment Calculate the current search point respectively perpendicular distance to the line segment This refers to the shortest distance. Among all these distances, the minimum value is taken as the current search point. Cost of controlling the distance to the corresponding initial inlay line .

[0138] To demonstrate the effectiveness of the technology proposed in this invention, we conducted the following experiments:

[0139] Qualitative evaluation

[0140] Figure 10 , 11 The proposed optimized mosaicking network is presented, along with the mosaicking results of a vector file for a specific region generated based on this network. We selected 10 local regions from the mosaicking results, see [link to documentation]. Figure 11 In regions 1 and 2, the mosaic lines follow the riverbank, avoiding islands in the river; in regions 3, 4, and 6, the mosaic lines successfully bypass ponds, farmland, and woodlands within villages with small grayscale differences; in region 8, the mosaic lines avoid cloud-covered areas; in regions 5, 10, 11, and 12, the mosaic lines follow the main roads, avoiding crossing buildings; and in regions 7 and 9, the mosaic lines follow areas near valleys, avoiding crossing ridges. These results demonstrate the strong robustness of the method.

[0141] In high-resolution urban image mosaicking, due to the dense buildings, complex road layouts, and the influence of side-view imaging angles, mosaic line extraction must comprehensively consider building edges, road grayscale differences, and texture features. (See...) Figure 12 In vertically captured images, the mosaic lines generated by the mosaic method of this invention can effectively avoid the main building in most positions and advance along the edge of a road or open space, see... Figure 12 (a)(1)(2)(3). However, for urban areas viewed from the side, when the mosaic line originates in a densely built-up area, sloping building structures may obstruct the road between buildings. Although the mosaic line attempts to follow the road, it may still pass through some buildings. In wide roads and open areas, the mosaic line effectively avoids buildings, see [reference needed]. Figure 12 (b)(2); and in areas with lower building heights, the mosaic lines also follow the road well (see...). Figure 12 (b)(3)). However, in Figure 12 In (c)(2)(3), due to the adaptive grayscale threshold weight setting adopted by this method, it is easily affected by the shadow of high-rise buildings, resulting in the mosaic line part crossing the building.

[0142] For mosaicking high-resolution urban imagery, when the imaging angle is close to vertical, this method can effectively utilize grayscale, edge, and texture information, allowing the mosaic lines to avoid buildings in most locations and advance along the boundaries of roads or open spaces. (See...) Figure 12(a) When a significant side-view imaging angle exists, especially in areas with dense high-rise buildings and severe sloping obstructions, buildings may block surface roads, causing the mosaic line to partially penetrate the building structure, see [reference needed]. Figure 12 (b)(2)(3) Figure 12 (c)(2)(3). Although the introduction of adaptive grayscale thresholds improves the recognition of roads and low-rise building areas, high-rise buildings and their shadows still interfere with the extraction of mosaic lines. Overall, this method performs better in vertically captured urban images; however, for images with large side-view angles, it is expected that the accuracy and stability under severe occlusion conditions can be further improved by incorporating building height or 3D model information.

[0143] Quantitative evaluation

[0144] To further evaluate the robustness of our method in these two typical regions, we compared images from high-density urban areas and low-density mountainous areas using existing large-scale image mosaicking methods and mainstream remote sensing image processing software, such as... Figure 13 , Figure 14 As shown. The comparison methods include the traditional A* algorithm, the SMP-DP algorithm, the ERDAS software proposed by Chai et al., and the Ortho-Vista software. First, we compare the running time of the different methods, as well as the number of protruding buildings intersecting in the same area, and the radial and geometrical consistency of the mosaic lines extracted using these methods.

[0145] exist Figure 13 and Figure 14 In this study, we compared the performance of mosaic lines from different methods in local areas. Both our method and traditional A methods can avoid most buildings and rural farmland, choosing to traverse roads or areas with small grayscale variations. In contrast, traditional A methods, due to the lack of structural similarity, perform weakly in areas with differences in illumination intensity between the left and right images. SMP-DP and ERDAS methods perform poorly in areas with discontinuous edge information and have limited path search directions, easily crossing buildings or rivers. Chai et al.'s method is more flexible, but still crosses prominent terrain.

[0146] Experiments show that our method traverses significantly fewer buildings than other methods, for example, by 14.29% compared to traditional A* and 20.75% compared to SMP-DP. In mountainous areas, our method avoids areas with significant terrain undulations and thin clouds in lake areas by using an adaptive obstacle threshold, while other methods tend to traverse areas with high grayscale variations.

[0147] In terms of operational efficiency, although SMP-DP has a shorter runtime, its results are inferior to our method in high-quality mosaic line extraction due to path direction limitations. Overall, our method performs better in avoiding buildings and complex terrain, and has higher mosaic quality. It is suitable for various scenarios such as urban plains, mountains, and lakes. Table 1 shows a comparison of the runtime and the number of obvious buildings traversed by each method in urban and rural areas.

[0148] Table 1 Quantitative Comparison

[0149]

[0150] Radiometric consistency assessment primarily evaluates color difference along the mosaic line. Smaller differences indicate a smoother color transition between adjacent images on either side of the mosaic line, resulting in less noticeable seams between images. Furthermore, we calculated the maximum, minimum, and average grayscale differences between adjacent images along the mosaic line path. The results are shown in the table below. The average grayscale difference between adjacent images along the mosaic line is calculated as follows:

[0151]

[0152] here Coordinates of the points along the mosaic path and Point The grayscale values ​​corresponding to the coordinates are shown in adjacent images. The average grayscale difference between the mosaic path points in the left and right images is calculated point-by-point using the formula described above. A smaller average value indicates better performance in extracting the mosaic path, while a larger average value indicates worse performance. The radiation consistency of each method in urban and rural areas is compared in Tables 2 and 3.

[0153] Table 2 Comparison of Radiation Consistency in Urban Areas

[0154]

[0155] Table 3 Comparison of Radiation Consistency in Mountain Village Areas

[0156]

[0157] Quantitative comparisons show that our method generates mosaic lines with the smallest maximum grayscale difference in densely built-up urban plains, with an average grayscale difference of 10.0164, which is 46.12% lower than traditional A, 59.97% lower than map cutting, 25.97% lower than Ortho-Vista, and 20.94% lower than ERDAS. In mountainous towns and lake areas, the average grayscale difference is 10.9016, which is 20.91% lower than A, 30.98% lower than map cutting, and 32.26% lower than ERDAS, comparable to SMP-DP and Ortho-Vista. Our method traverses the fewest protruding buildings, reducing misalignment caused by feature variations or errors. The optimized mosaic lines generate superior mosaic images, meeting the needs of large-scale digital orthophoto stitching.

[0158] Geometric consistency evaluation employs feature matching. A schematic diagram of the geometric consistency evaluation method can be found here. Figure 15 Due to limitations in exterior orientation elements and DEM accuracy, geometric misalignment may still exist in the corrected orthophoto. To evaluate the geometric consistency of mosaic lines extracted by different methods, this study divided the overlapping region into 500×500 pixel blocks, used GLS-MIFT technology to match feature points, eliminated mismatches, and integrated the matching points of all blocks. A 10-pixel wide buffer was constructed around the mosaic lines, and feature matching points were extracted from the buffers of the two images used for mosaicking, referred to here as the left image or the right image, respectively. Figure 15 (a)(b). Connecting the matching points in the buffer represents the common feature points in the buffer; these common feature points are the matching results. See [link to relevant documentation]. Figure 15 (c) The GLS-MIFT algorithm identifies more feature points than the traditional SIFT algorithm, ensuring a sufficient number of feature points within the buffer. Based on this, the degree of geometric misalignment within the mosaic line buffer is quantitatively measured. The absolute difference in the x and y directions of all matching feature points and the displacement distance D are calculated. x and y represents the difference between the x and y directions, and D can be represented as:

[0159]

[0160] Next, we calculated the average values ​​of the row and column coordinate differences and displacement distances within the overlapping regions. For statistical fairness, at least 500 pairs of feature matching points were retained in the buffer constructed from the mosaic lines generated by each method, after excluding feature matching points with displacement distances greater than 4. Average values , and Calculate using the following formula:

[0161]

[0162]

[0163]

[0164] In these formulas and Let x and y represent the differences between the x and y coordinates of the i-th feature matching point in the left and right image buffers, respectively. This represents the displacement distance of the i-th feature matching point within the buffer, and n represents the number of feature matching points within the buffer. By comparison... , and These values ​​can be used to assess the geometric quality of image mosaicking. The smaller these values, the better the geometric quality, and vice versa.

[0165] Table 4 Comparison of Geometric Consistency in Urban Areas

[0166]

[0167] Table 5. Comparison of Geometric Consistency in Mountain Village Areas

[0168]

[0169] Tables 4 and 5 show that the method of this invention performs well, followed by graph cutting, traditional A*, SMP-DP, and ERDAS. The method of this invention performs feature detection based on SSIM, ensuring geometric consistency and performing well even in sparsely textured mountainous areas. Chai et al.'s method improves the initial mosaic line through pixel-level search, resulting in more stable performance. SMP-DP and ERDAS, limited by edge information and grayscale difference weight settings, perform poorly in terms of geometric consistency.

[0170] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An orthophoto mosaicking method using the adaptive cost A* algorithm, characterized in that... Includes the following steps: Step 1: Obtain the initial tessellation lines between each pair of orthophoto images to be tessellated, so as to construct an initial tessellation line network based on the polygon skeleton; Step 2: Optimize the initial mosaicking line network using an improved A* algorithm; An adaptive threshold algorithm is used to determine the passable area and construct a cost function to obtain an improved A* algorithm. Then, the improved A* algorithm is used to optimize the initial mosaic lines in the initial mosaic line network one by one to obtain the optimized mosaic lines, thereby completing the mosaicking between orthophoto images. The global threshold is the average gray value of the overlapping area of ​​adjacent images. The local threshold is determined by the average gray value of N*N blocks in the left and right images where the current search point is located. The adaptive threshold algorithm adaptively selects either the global threshold or the local threshold as the adaptive threshold and compares it with the maximum gray value difference cost of the current search point to determine whether the current search point is a passable area. When using the improved A* algorithm for searching, the adaptive threshold algorithm is first used to determine whether the current search point is a passable area. If not, the search proceeds to the next pixel. If so, then search the pixels in the eight-directional neighborhood around the current search point. If any pixel satisfying formulas (5) and (6) is found, then the grayscale cost is obtained. , （5） （6） Then, using the following formula, the cost function F is calculated, and subsequent search steps are performed. Where D represents distance control cost, G represents the distance from the starting point to the current search point, H represents the distance from the current search point to the destination, w represents the estimated path cost weight, Grad represents the grayscale cost, c represents the grayscale conversion coefficient used to balance the weight settings, and SSIM represents the structural similarity index. Indicates starting from the current search point The Euclidean distance to the starting point, where α represents a positive integer greater than 1, and T represents the adaptive threshold. This represents the cost of the maximum grayscale difference. This represents the cost of the intermediate grayscale difference; If no pixel satisfying formulas (5) and (6) is found, the search continues until the search is completed and the optimized mosaic line is obtained. When using the improved A* algorithm for the search, the start and end points of the initial tessellation line are adjusted using the following method. Based on the current search point Adaptive threshold Intermediate grayscale difference cost The following conditions must be met for the determination: like Then determine the current search point. If it is a non-obstacle point, then it is an obstacle point; otherwise, it is an obstacle point. With the center as The process involves traversing pixels with a radius of dN. If all pixels are obstacles, then the current search point is... Passable depth ; If both conditions are met and This indicates the current search point. These are non-obstacle points and require no adjustment. Otherwise, perform a depth-first traversal of the surrounding neighborhood until a satisfactory result is found. and The pixels are used to replace the current search point. .

2. The orthophoto mosaicking method using the adaptive cost A* algorithm according to claim 1, characterized in that: When using the adaptive threshold algorithm to determine passable areas, the adaptive threshold T is first calculated using the following formula. （1） （2） in, The global threshold represents the average grayscale value of the overlapping region between adjacent images. These represent the current search points. The average grayscale value of the N*N blocks in the adjacent left and right images is calculated using formula (2). Indicates the current search point grayscale value; from The local threshold is selected based on the minimum grayscale average value. If the local threshold is less than the global threshold, the local threshold is used as the adaptive threshold T; otherwise, the global threshold is used as the adaptive threshold T. Then, using the following formula, calculate the maximum grayscale difference cost corresponding to the current search point. （3） （4） in, These represent the current search points in adjacent left and right images, respectively. grayscale value, These represent the current search points in adjacent left and right images, respectively. The maximum gray-level gradient in the eight directional neighborhoods: top left, top, top right, right, bottom right, bottom, bottom left, and left. Finally, if the current search point Corresponding maximum grayscale difference cost If the value is greater than the adaptive threshold T, the pixel is determined to be a non-passable area; otherwise, the pixel is determined to be a passable area.

3. The orthophoto mosaicking method using the adaptive cost A* algorithm according to claim 1, characterized in that: Calculate the distance control cost D using the following formula. （9） in, This represents any line segment between two adjacent points on the initial tessellation line. Indicates the current search point to line segment The vertical distance.

4. The orthophoto mosaicking method using the adaptive cost A* algorithm according to claim 1, characterized in that: When constructing the initial mosaicking network, the orthophoto image to be mosaicked first needs to be preprocessed, including geometric correction and downsampling.

Images