DEM image slicing method based on multi-gpu computing mode
By employing a multi-GPU parallel computing mode and CUDA stream overlay technology, the problem of excessively long processing time for traditional CPUs in DEM image processing has been solved, achieving efficient DEM image tiling and meeting the real-time processing needs of massive amounts of data.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA RAILWAY DESIGN GRP CO LTD
- Filing Date
- 2025-12-31
- Publication Date
- 2026-05-29
AI Technical Summary
Traditional CPU-based serial data processing methods struggle to meet real-time requirements when processing massive amounts of DEM image data, and the insufficient computing power of a single GPU leads to low efficiency in processing high-resolution DEM images.
By employing a multi-GPU parallel computing mode, combined with CUDA stream overlay technology and the Ring-Allreduce communication algorithm, and through asynchronous transmission and parallel computing optimization, efficient slicing processing of DEM images is achieved.
It significantly improves the processing efficiency of DEM image slices, solves the problems of excessive time consumption and insufficient video memory in traditional methods, and realizes real-time processing of large-scale, high-resolution DEM data.
Smart Images
Figure CN122115189A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of geospatial data processing technology, specifically relating to a DEM image tiling method based on a multi-GPU computing mode. Background Technology
[0002] With the rapid development of geospatial data acquisition technology, the scale of Digital Elevation Model (DEM) data has exploded. Traditional CPU-based serial data processing methods struggle to meet the demands for real-time data processing and efficient analysis when faced with massive amounts of data. Specifically, CPUs have a relatively small number of computing cores, and each core has limited computing power, making it difficult to fully utilize the advantages of multi-core architectures. When processing high-resolution DEM image data, CPUs need to perform calculations pixel by pixel, resulting in excessively long processing times and failing to meet real-time requirements.
[0003] In recent years, with the rapid development of GPU technology, its application in the field of parallel computing has gradually attracted attention. GPUs have a large number of computing cores, enabling them to efficiently handle large-scale parallel computing tasks. However, the computing power and memory capacity of a single GPU are still limited, and for larger-scale DEM image data processing tasks, the performance bottleneck of a single GPU is becoming increasingly apparent. Therefore, how to fully utilize the parallel computing capabilities of multiple GPUs to further improve the efficiency of DEM image tiling processing has become an urgent problem to be solved. Summary of the Invention
[0004] To address the shortcomings of existing technologies, this invention provides a DEM image tiling method based on a multi-GPU computing mode that can further improve processing efficiency and meet the needs of real-time processing of massive amounts of data.
[0005] Therefore, the present invention adopts the following technical solution:
[0006] A DEM image tiling method based on multi-GPU computing mode includes the following steps:
[0007] S1, Obtain the features of the DEM image file and form a feature table. The polygon is obtained from the input polygon point set. ;
[0008] S2, Calculate the polygon The intersection points of each edge of the polygon with the bounding rectangle of each DEM image are merged to form a complete polygon. Intersection information with DEM imagery, and will be fully integrated into polygons. DEM image data within the range is added to the intersection information;
[0009] S3 evenly distributes the DEM image across the GPUs, starts an independent stream for each GPU, allocates video memory, creates a mask and transmits data asynchronously, thus achieving the overlap of computation and communication.
[0010] S4, starts a kernel function on each GPU to perform slicing operations on the allocated DEM image to obtain the DEM image within the polygon range, and merges the DEM image slices.
[0011] S5 copies the merged DEM image slice processing results back to CPU memory and outputs the results.
[0012] In step S1 above, the features of the DEM image file are obtained to form a feature table. The specific operation is as follows: sequentially traverse the disk directories on the host side. For each DEM image, its metadata is read to obtain the coordinates of the lower left corner of the DEM image. upper right corner coordinates Pixel resolution and number of rows and columns ; Calculate the bounding rectangle of the DEM image The four vertices and the center point The coordinates; the storage path and Composition of feature tuples Ultimately, a feature table is formed. .
[0013] In step S1 above, a polygon is obtained. The specific operation is as follows: Read in the unordered point set of the polygon. , To find the number of vertices of the polygon, first find its centroid. Then calculate the polar angle. ; and according to Sort the polygons in ascending order; if self-intersections occur after sorting, automatically switch to the Graham scan algorithm to generate the convex hull, ensuring the final polygon... There are no self-intersections and the vertex order is counterclockwise.
[0014] The specific steps of step S2 above are as follows:
[0015] S21, the polygon Disassembled into sides The set of DEM feature tables obtained from S1 Extract all bounding rectangles of the DEM. polygon Each edge With all DEM bounding rectangles intersection Solve using linear algebra, the formula is as follows:
[0016] ,
[0017] in, It is a scalar parameter in the interval [0,1].
[0018] All calculated intersection points Corresponding DEM image index and polygons Each edge index Record them to form an intersection table. ;
[0019] S22, merge the intersection results of all edges, indexed by DEM image. Aggregate intersection table Each DEM image and polygon are obtained. Intersection information description This forms complete intersection information between the polygon and the DEM image;
[0020] S23, Add DEM image data that is entirely within the polygon area; if the bounding rectangle of a DEM image... Completely located within a polygon Inside, then The four vertices are directly added to the intersection information description. .
[0021] The specific steps of step S3 above are as follows:
[0022] S31, will The DEM image is uniformly mapped to the hash-partitioned image using a hash-partitioning strategy. On a GPU, Responsible for index set , making ;
[0023] S32, for Creating non-blocking CUDA streams This ensures that subsequent data transfer and computation overlap; memory blocks are pre-allocated based on the maximum slice size. Size, of which The maximum number of rows in all DEMs to be processed is The maximum number of columns is The size of the allocated video memory is ;
[0024] S33, Creation and the DEM image Binary polygonal mask with perfectly aligned grid sizes ,when The grid cells fall into the polygon Inside, Otherwise, it is 0; where, and represents the x and y coordinates of the pixel;
[0025] S34, employs asynchronous transmission technology to combine DEM images and polygon masks. Parallel transmission to each GPU memory; asynchronous transmission is a GPU computing optimization method, where the CPU transmits data to the GPU, and the transmission process is synchronized with the calculation of the transmitted data in S4.
[0026] The specific steps of step S4 above are as follows:
[0027] S41, slice each DEM image; The kernel function `clip_kernel` is launched, and each thread processes DEM image pixels. First, calculate the geographic coordinates according to the affine transformation of the target. Then, the source pixel index is calculated in reverse. If the source pixel index If the image falls within the spatial range of the DEM to be tiled, then the DEM image to be tiled will be indexed at the source pixel index. elevation value at the location With polygonal mask Multiplication yields a temporary value After traversing all pixels of a single DEM, the single DEM slice is completed.
[0028] S42, merge slices from a single DEM image; immediately after slicing, use the atomic operation cuda.atomic.max to merge the slices. With video memory blocks Target pixel position Perform atomic maximum value reduction, and then loop through all threads. The internal storage content is The merged result of all DEM image slices within the polygon range, achieving integrated slicing and merging;
[0029] S43 synchronizes the CUDA execution flow of all GPUs and merges the slice processing results of each GPU;
[0030] Call cudaStreamSynchronize( The process waits for all GPUs' CUDA execution flows to complete, and then merges local results between GPUs using the Ring-AllReduce ring topology communication algorithm: In the Ring-AllReduce ring communication topology, each GPU sends only a portion of its data block to the next adjacent GPU, while simultaneously receiving another portion of data from the previous GPU; the GPUs synchronously exchange and merge data, and perform element-level maximum operations to preserve the maximum value in overlapping areas; after k-1 iterations, each GPU can obtain the global DEM image tile matrix. .
[0031] The specific operation of step S5 above is as follows: Select any one Use cudaMemcpyAsync to slice the global DEM image matrix. The data is asynchronously copied from GPU memory back to CPU memory and written to a standard GeoTIFF file for output.
[0032] Compared with the prior art, the present invention has the following beneficial effects:
[0033] 1. This invention is the first to propose using a multi-GPU parallel architecture in CUDA mode for fragmented DEM tiling, which solves the problems of excessive time consumption for traditional CPU pixel-by-pixel processing and insufficient single GPU video memory, which prevents real-time high-resolution terrain tiling.
[0034] 2. This invention employs asynchronous transmission and CUDA stream overlap technology to enable parallel data transmission on the PCIe bus and GPU kernel function computation. By allocating data copy operations and computation tasks to different CUDA streams, the CPU can initiate the transmission of the next piece of data while the GPU is performing computation, avoiding GPU idleness due to waiting for data transmission. This improves PCIe bandwidth and GPU computing core utilization, achieving higher overall throughput efficiency.
[0035] 3. This invention utilizes the Ring-Allreduce architecture to enable computing nodes to scale synchronously with the amount of data: For larger-scale, higher-resolution DEM data, tasks can be split into finer-grained parts and evenly distributed to more GPUs; at the same time, the ring-Allreduce's circular communication mechanism, which only communicates with adjacent nodes, ensures that data transmission time increases slowly with the number of GPUs, without offsetting the computational speedup. Therefore, the speedup ratio can increase linearly with the scale of GPUs, effectively alleviating the performance bottleneck of large-scale DEM processing.
[0036] 4. This invention effectively alleviates the performance bottleneck of large-scale DEM processing by using the Ring-Allreduce communication algorithm:
[0037] For large-scale, high-resolution DEM data, tasks are split into smaller, finer-grained parts and evenly distributed to more GPUs for processing. Meanwhile, the Ring-Allreduce communication algorithm, which communicates only with adjacent nodes, reduces the data transfer time consumption that increases with the number of GPUs, without offsetting the computational speedup from the increase in the number of GPUs. The speedup ratio can increase linearly with the scale of GPUs. Attached Figure Description
[0038] Figure 1 This is a flowchart of a DEM image tiling method based on multi-GPU computing mode in an embodiment of the present invention;
[0039] Figure 2 The diagram shows the process of polygon tiling of DEM image based on multi-GPU computing mode, where (a) is the original strip-shaped discrete distribution of DEM image block, (b) is the cropped polygon that defines the DEM tiling range, and (c) is the complete and seamless DEM polygon tiling output after multi-GPU parallel tiling and merging.
[0040] Figure 3 This is a comparison chart showing the efficiency of multi-GPU parallel computing and traditional serial computing under different data scales. Detailed Implementation
[0041] The technical solution of the invention will be clearly and completely described below with reference to the accompanying drawings and embodiments. Obviously, the following embodiments are only some embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0042] Example
[0043] This invention provides a DEM image tiling method based on a multi-GPU computing mode, such as... Figure 1 As shown, the specific steps are as follows:
[0044] S1: Obtain the file features of the DEM image and create a polygon list based on the input polygon point set. Details are as follows:
[0045] S11, Scan the input DEM dataset and extract the spatial feature information of each DEM image file.
[0046] On the host side, sequentially traverse the disk directories. Zhang DEM image, Figure 2 (a) is shown as The DEM image is composed of a strip-shaped discrete distribution of DEM images. Its metadata is read (e.g., using rasterio.open) to obtain the coordinates of the lower left corner of the DEM image. upper right corner coordinates Pixel resolution and number of rows and columns Calculate the bounding rectangle of the DEM image. The four vertices and the center point The coordinates; the storage path and Composition of feature tuples Ultimately, a feature table is formed. .
[0047] S12, for the input polygon point set, arrange them sequentially into a planar polygon, such as... Figure 2 As shown in the red box in (b), a list of polygons is formed. The specific operation is as follows: read in the unordered set of points of the polygons. , To find the number of vertices of the polygon, first find its centroid. Then calculate the polar angle. ; and according to Sort the polygons in ascending order; if self-intersections occur after sorting, automatically switch to the Graham scan algorithm to generate the convex hull, ensuring the final polygon... There are no self-intersections and the vertex order is counterclockwise.
[0048] S2, calculate the intersection points of each side of the polygon with the bounding rectangle of each DEM image, merge them to form complete intersection information between the polygon and the DEM image, and add the DEM image data that is completely within the polygon's range. The specific steps are as follows:
[0049] S21, the polygon Disassembled into sides The set of DEM feature tables obtained from S11 Extract all bounding rectangles of the DEM. polygon Each edge With all DEM bounding rectangles intersection The solution can be obtained using linear algebra, as shown in the following formula:
[0050]
[0051] in, It is a scalar parameter in the interval [0,1].
[0052] All intersection points obtained from the calculation Corresponding DEM image index and polygons Each edge index Record them to form an intersection table. .
[0053] S22, merge the intersection results of all edges, according to the DEM index. Aggregate intersection table Each DEM image and polygon are obtained. Intersection information description This forms the intersection information between the complete polygon and the DEM image.
[0054] S23, Add DEM image data that is entirely within the polygon area; if the bounding rectangle of a DEM image... Completely located within a polygon Inside, then The four vertices are directly added to the intersection information description. This avoids redundant calculations later.
[0055] S3 distributes the DEM evenly across the GPUs, starts an independent stream for each GPU, allocates memory size, creates a mask, and asynchronously transmits data, achieving overlap between computation and communication. The specific steps are as follows:
[0056] S31, will The DEM image is uniformly mapped to the hash-partitioned image using a hash-partitioning strategy. On a GPU, Responsible for index set , making .
[0057] S32, which starts an independent CUDA execution stream for each GPU; that is... Creating non-blocking CUDA streams This ensures that subsequent data transfer and computation can be performed concurrently. Memory blocks are pre-allocated based on the maximum slice size. Size, of which The maximum number of rows in all DEMs to be processed is The maximum number of columns is The size of the allocated video memory is .
[0058] S33, create a polygonal mask, that is, create a mask that is similar to the first polygonal mask. DEM image Binary polygonal mask with perfectly aligned grid sizes ,when The grid cells fall into the polygon Inside, Otherwise, it is 0. Wherein, and represents the x and y coordinates of the pixel.
[0059] S34, employs asynchronous transmission technology to combine DEM images and polygon masks. Parallel transfer to each GPU's video memory. Asynchronous transfer is a GPU computing optimization method, where the CPU transfers data to the GPU, and the transfer process is synchronized with the calculation of the transferred data in S4.
[0060] S4 starts a kernel function on each GPU to perform slicing operations on the allocated DEM image to obtain the DEM image within the polygon range, and merges the DEM image slices.
[0061] S41, slice each DEM image; The kernel function `clip_kernel` is launched, and each thread processes DEM image pixels. First, calculate the geographic coordinates according to the affine transformation of the target. Then, the source pixel index is calculated in reverse. If the source pixel index If the image falls within the source boundary (the spatial range of the DEM image to be tiled), then the source elevation (the DEM image at the source pixel index) will be used. (elevation value at the location) With polygonal mask Multiplication yields a temporary value After traversing all pixels of a single DEM, the single DEM slice is completed.
[0062] S42, merge slices from a single DEM image; immediately after slicing, use the atomic operation cuda.atomic.max to merge the slices. With video memory blocks Target pixel position Perform atomic maximum value reduction, and then loop through all threads. The stored content is The results of merging all DEM image slices within a polygonal range are integrated to achieve unified slicing and merging.
[0063] S43 synchronizes the CUDA execution flow of all GPUs and merges the slice processing results of each GPU.
[0064] Call cudaStreamSynchronize( After all GPUs' CUDA execution flows complete, local results are merged between GPUs using the Ring-AllReduce algorithm: In the Ring-AllReduce topology, each GPU sends a portion of its data block to the next adjacent GPU while receiving another portion from the previous GPU; GPUs synchronously exchange and merge data, and perform element-level maximum operations to preserve the maximum value in overlapping regions; after k-1 iterations (k being the total number of GPUs), each GPU obtains the global DEM slice matrix. ,like Figure 2 As shown in (c).
[0065] S5 copies the merged DEM slice processing results back to CPU memory and outputs the results.
[0066] Choose any one Use cudaMemcpyAsync to slice the global DEM image matrix. The data is asynchronously copied from GPU memory back to CPU memory and written to a standard GeoTIFF file for output.
[0067] like Figure 3 As shown, at a resolution of 30m, the CPU takes 68.98s to process DEM image tiles, while the GPU takes 11.28s, representing a speedup of up to 6 times. This demonstrates that the method of this invention significantly improves processing efficiency. When the image resolution increases from 200m to 30m and the data volume increases 40 times, the speedup ratio increases from 2.34× to 6.12×. As can be seen from the figure, the speedup ratio increase in the high-resolution range is significantly higher than that in the low-resolution range, verifying that the larger the data volume, the more significant the performance gain brought by GPU expansion. This proves that the method of this invention can effectively alleviate the performance bottleneck in large-scale DEM processing.
[0068] In summary, this invention significantly improves the processing efficiency of DEM image tiling by using an efficient parallel tiling algorithm and a multi-GPU collaborative processing mechanism, while also optimizing the data distribution and transmission strategy through asynchronous transmission and CUDA streaming, thus meeting the needs for real-time processing of massive DEM data.
Claims
1. A DEM image tiling method based on multi-GPU computing mode, characterized in that, Includes the following steps: S1, Obtain the features of the DEM image file and form a feature table. The polygon is obtained from the input polygon point set. ; S2, Calculate the polygon The intersection points of each edge of the polygon with the bounding rectangle of each DEM image are merged to form a complete polygon. Intersection information with DEM imagery, and will be fully integrated into polygons. DEM image data within the range is added to the intersection information; S3 evenly distributes the DEM image across the GPUs, starts an independent stream for each GPU, allocates video memory, creates a mask and transmits data asynchronously, thus achieving the overlap of computation and communication. S4, starts a kernel function on each GPU to perform slicing operations on the allocated DEM image to obtain the DEM image within the polygon range, and merges the DEM image slices. S5 copies the merged DEM image slice processing results back to CPU memory and outputs the results.
2. The DEM image slicing method according to claim 1, characterized in that: In S1, features of the DEM image file are obtained to form a feature table. The specific operation is as follows: sequentially traverse the disk directories on the host side. For each DEM image, its metadata is read to obtain the coordinates of the lower left corner of the DEM image. upper right corner coordinates Pixel resolution and number of rows and columns ; Calculate the bounding rectangle of the DEM image The four vertices and the center point The coordinates; the storage path and Composition of feature tuples Ultimately, a feature table is formed. .
3. The DEM image slicing method according to claim 1, characterized in that: In S1, we obtain the polygon. The specific operation is as follows: Read in the unordered point set of the polygon. , To find the number of vertices of the polygon, first find its centroid. Then calculate the polar angle. ; and according to Sort the polygons in ascending order; if self-intersections occur after sorting, automatically switch to the Graham scan algorithm to generate the convex hull, ensuring the final polygon... There are no self-intersections and the vertex order is counterclockwise.
4. The DEM image slicing method according to claim 3, characterized in that, The specific steps of S2 are as follows: S21, the polygon Disassembled into sides The set of DEM feature tables obtained from S1 Extract all bounding rectangles of the DEM. polygon Each edge With all DEM bounding rectangles intersection Solve using linear algebra, the formula is as follows: ， in, It is a scalar parameter in the interval [0,1]. All calculated intersection points Corresponding DEM image index and polygons Each edge index Record them to form an intersection table. ; S22, merge the intersection results of all edges, indexed by DEM image. Aggregate intersection table Each DEM image and polygon are obtained. Intersection information description This forms complete intersection information between the polygon and the DEM image; S23, Add DEM image data that is entirely within the polygon area; if the bounding rectangle of a DEM image... Completely located within a polygon Inside, then The four vertices are directly added to the intersection information description. .
5. The DEM image slicing method according to claim 1, characterized in that, The specific steps for S3 are as follows: S31, will The DEM image is uniformly mapped to the hash-partitioned image using a hash-partitioning strategy. On a GPU, Responsible for index set , making ; S32, for Creating non-blocking CUDA streams This ensures that subsequent data transfer and computation overlap; memory blocks are pre-allocated based on the maximum slice size. Size, of which The maximum number of rows in all DEMs to be processed is The maximum number of columns is The size of the allocated video memory is ; S33, Creation and the DEM image Binary polygonal mask with perfectly aligned grid sizes ,when The grid cells fall into the polygon Inside, Otherwise, it is 0; where, and represents the x and y coordinates of the pixel; S34, employs asynchronous transmission technology to combine DEM images and polygon masks. Parallel transmission to each GPU memory; asynchronous transmission is a GPU computing optimization method, where the CPU transmits data to the GPU, and the transmission process is synchronized with the calculation of the transmitted data in S4.
6. The DEM image slicing method according to claim 5, characterized in that, The specific steps for S4 are as follows: S41, slice each DEM image; The kernel function `clip_kernel` is launched, and each thread processes DEM image pixels. First, calculate the geographic coordinates according to the affine transformation of the target. ; Then calculate the source pixel index. If the source pixel index If the image falls within the spatial range of the DEM to be tiled, then the DEM image to be tiled will be indexed at the source pixel index. elevation value at the location With polygonal mask Multiplication yields a temporary value After traversing all pixels of a single DEM, the single DEM slice is completed. S42, merge single DEM image slices; Immediately after slicing, use the atomic operation cuda.atomic.max to... With video memory blocks Target pixel position Perform atomic maximum value reduction, and then loop through all threads. The internal storage content is The merged result of all DEM image slices within the polygon range, achieving integrated slicing and merging; S43 synchronizes the CUDA execution flow of all GPUs and merges the slice processing results of each GPU; Call cudaStreamSynchronize( The process waits for all GPUs' CUDA execution flows to complete, and then merges local results between GPUs using the Ring-AllReduce ring topology communication algorithm: In the Ring-AllReduce ring communication topology, each GPU sends only a portion of its data block to the next adjacent GPU, while simultaneously receiving another portion of data from the previous GPU; the GPUs synchronously exchange and merge data, and perform element-level maximum operations to preserve the maximum value in overlapping areas; after k-1 iterations, each GPU can obtain the global DEM image tile matrix. .
7. The DEM image slicing method according to claim 6, characterized in that, The specific operation of S5 is as follows: Select any card Use cudaMemcpyAsync to slice the global DEM image matrix. The data is asynchronously copied from GPU memory back to CPU memory and written to a standard GeoTIFF file for output.