A method, device and medium for stitching multispectral images of unmanned aerial vehicle remote sensing

Abstract

The embodiment of the present invention discloses a method, device and medium for stitching multispectral images of remote sensing by unmanned aerial vehicles. The method comprises: acquiring N frames of single-band images and performing band registration correction to obtain N frames of multispectral images, wherein each frame of single-band image contains source images of M bands; performing band fusion on the N multispectral images to obtain N grayscale images, and stitching the N grayscale images to obtain N local homography matrix sets; and obtaining the local homography matrix set D by the nth local homography matrix set. n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band spliced ​​images, and the M single-band spliced ​​images are fused to obtain a preliminary spliced ​​image; the line segment structure existing in the preliminary spliced ​​image is inspected and corrected to obtain the final splicing result.

Description

Technical Field

[0001] The present invention relates to the field of image processing, and more specifically, to a method, device and medium for stitching multispectral images of remote sensing by unmanned aerial vehicles. Background Art

[0002] With the continuous development of spectral imaging technology, drones equipped with lighter spectral imagers have been widely used in geographic remote sensing, agricultural remote sensing and other fields. Since drones fly at a low altitude, the image coverage area collected each time is limited. In order to obtain more comprehensive information about the entire scene, multiple images need to be stitched together to obtain a panoramic image with a larger viewing angle. However, the existing stitching methods are more targeted at full-color images or RGB images, and there are fewer stitching methods for multispectral images; and there will be distortion and deformation during the stitching process. Due to the large number of stitched images and the uneven image quality, most current algorithms are difficult to obtain efficient and correct stitching results. Summary of the invention

[0003] The object of the present invention is to provide a method, device and medium for stitching multispectral images of unmanned aerial vehicle remote sensing, so as to solve at least one of the problems existing in the prior art.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] The first aspect of the present invention provides a method for stitching multispectral images of unmanned aerial vehicle remote sensing, comprising:

[0006] S10: Acquire N frames of single-band images and perform band registration correction to obtain N frames of multispectral images, wherein each frame of single-band image contains source images of M bands, wherein:

[0007] N and M are positive integers;

[0008] S20: performing band fusion on the N multispectral images to obtain N grayscale images, and performing splicing based on the N grayscale images to obtain N local homography matrix sets;

[0009] S30: Through the nth local homography matrix set D n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band spliced ​​images, and the M single-band spliced ​​images are fused to obtain a preliminary spliced ​​image, wherein,

[0010] n=1, 2,...,N; m=1, 2,...,M;

[0011] S40: Checking and correcting the line segment structure in the preliminary stitched image to obtain a final stitching result.

[0012] Further, the S10 includes:

[0013] S101: acquiring a video along a specific route through a multispectral camera of an unmanned aerial vehicle, wherein the video includes N frames of single-band images, and each frame of the single-band image includes source images of M bands;

[0014] S102: performing band registration correction on the N frames of single-band images respectively by using the SIFT algorithm: performing convolution with the current single-band image by using Gaussian kernel functions of different sizes, reducing the image after the convolution at different levels to obtain a scale space, and performing extreme value detection on the source matrix and the mesh matrix in the scale space to obtain feature points, while removing edge response points in the feature point concentration, and assigning direction values ​​to key points;

[0015] Generate a feature point descriptor according to the direction value, and match the feature points to obtain matching points;

[0016] The matching points are used to calculate the transformation matrix to perform image transformation and difference on the source images of the M bands of the current single-band image, obtain the multispectral image after band registration, and update the source images of the M bands.

[0017] Furthermore, the feature descriptor includes the scale, position and direction of the feature point.

[0018] Further, the S20 includes:

[0019] S201: performing band fusion on the N multispectral images to obtain N grayscale images;

[0020] S202: performing single-band registration on the N grayscale images respectively;

[0021] S203: Adopting an approximate projection transformation algorithm to stitch the N grayscale images to obtain corresponding N local transformation matrix sets.

[0022] Further, the S203 includes:

[0023] S2031: extracting and matching features of the N grayscale images respectively to obtain a pre-transformation matrix;

[0024] S2032: Divide the grayscale image into small grids evenly. The moving line linear transformation algorithm calculates the weighted square sum of the distances between the matching points in the current grid and the center of the grid to obtain the local homography transformation matrix set D. n ,in,

[0025] K is the total number of grids in the current grayscale image;

[0026] S2033: transforming the grayscale image using the local homography transformation matrix set, and fusing and splicing the transformed image with the reference image using a weighted average fusion method.

[0027] Furthermore, the local homography matrix of the kth grid of the nth grayscale image is The calculation formula is shown in formula (1). The local transformation matrix is ​​estimated by solving the eigenvector corresponding to the minimum eigenvalue Value:

[0028]

[0029] Among them, a i is the linear mapping function between matching points, l is the number of matching points; d is the distance between the matching point and the center of the grid, is the weight corresponding to the feature point, W k is the diagonal matrix corresponding to the weight, and A is the mapping matrix form corresponding to the mapping function, where:

[0030]

[0031]

[0032] Among them, x k and k is the coordinate of the center point of the kth grid, and is the coordinate of the i-th matching point in the k-th grid, σ is the Gaussian scaling factor, which ranges between (0,1).

[0033] Further, the S40 includes:

[0034] S401: Detect and extract line segment structures, extract line support regions in the preliminary spliced ​​image, and obtain line segment features through the line support regions, wherein:

[0035] The line support region is composed of pixel regions where line segments exist;

[0036] S402: Perform line segment correction: Use the least squares method to fit a line segment in the rectangular area, use the intersection of the line segments as feature points, sort their positional relationships with the center point of the rectangular area, calculate the perspective transformation matrix based on the correspondence between the transformed feature points and the original feature points, map all points in the area to be corrected, and complete the line segment correction.

[0037] Further, the S401 includes:

[0038] S4011: Use the four pixels to the lower right of each pixel to calculate the gradient corresponding to the pixel to obtain the position of the edge of the line segment. The gradient h of the pixel along the x direction and the y direction x (x,y) and h y The calculation formulas for (x, y) are shown in equations (4) and (5):

[0039]

[0040]

[0041] Where h(x, y) represents the grayscale value corresponding to the pixel with coordinates (x, y);

[0042] The calculation formula of gradient size H(x,y) and direction θ(x,y) is shown in formula (6):

[0043]

[0044] Among them, h x (x,y) and h y (x, y) are the gradients of the pixel along the x and y directions respectively;

[0045] S4012: Perform pseudo gradient sorting: divide the gradient value into 1024 equal parts from large to small, and take out a seed pixel from each segment for sorting;

[0046] Calculate the direction of the candidate rectangular area formed by the pixel points, and determine whether the pixel point belongs to the rectangular area by comparing the direction of the candidate pixel point and the candidate rectangular area, thereby completing the rectangular area division;

[0047] S4013: performing approximate calculation of rectangular areas: calculating the approximate rectangular areas corresponding to the line segments, the size of the rectangular areas covering the entire line support area, and each rectangular area represents a candidate line segment;

[0048] The center point c of the rectangular area x The calculation formula is shown in formula (7):

[0049]

[0050] Where (c x ,c y ) represents the coordinates of the center point, Region represents the current region, G(j) represents the gradient value of pixel j in the region, subscript j is used to traverse all pixels in the region, x(j) and y(j) represent the horizontal and vertical coordinates of pixel j respectively;

[0051] S4014: Improve the rectangular area and select the correct line support area: Calculate the number of false positives NFA value in the rectangular area. If the calculated result is less than a preset threshold, the current rectangular area is considered to be the correct line support area. The calculation formula of the NFA value is shown in formula (8):

[0052] NFA(r,t)=(PQ) 2.5 γB(f(r),k(r,t),p) (8)

[0053] Among them, P and Q are the image resolution sizes, B() represents the binomial distribution, f(r) is the number of pixels in the rectangular area r, k(r,t) represents the number of points within the class in the rectangular area r, p is the probability of gradient angle alignment, and γ is the possible number of different values ​​of p.

[0054] The second aspect of the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the method for stitching multispectral images of remote sensing of unmanned aerial vehicles provided in the first aspect of the present invention is implemented.

[0055] The third aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for stitching multispectral images of remote sensing of unmanned aerial vehicles provided in the first aspect of the present invention.

[0056] The beneficial effects of the present invention are as follows:

[0057] The solution provided in the present application first performs band registration correction on the grayscale images of each band to ensure the consistency between the bands, then obtains a set of local homography matrices by registering the fused single-band images to improve the efficiency of the algorithm, and finally performs line segment detection and correction on the spliced ​​image based on the line segment structure characteristics in the image to reduce the distortion and deformation that may be caused by splicing and improve the accuracy and precision of the splicing results. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] The specific implementation modes of the present invention will be further described in detail below in conjunction with the accompanying drawings.

[0059] Figure 1 A flowchart of a method for stitching multispectral images of remote sensing by an unmanned aerial vehicle provided by an embodiment of the present invention is shown.

[0060] Figure 2 A schematic diagram of a method for stitching multispectral images of UAV remote sensing provided by another embodiment of the present invention.

[0061] Figure 3 An embodiment of the present invention provides images of various bands before band registration.

[0062] Figure 4 An embodiment of the present invention provides images of each band after band registration.

[0063] Figure 5 A schematic diagram showing the division of a reference image provided by an embodiment of the present invention.

[0064] Figure 6 A schematic diagram showing the structure of a computer device for implementing the method provided by an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0065] In order to explain the present invention more clearly, the present invention is further described below in conjunction with embodiments and drawings. Similar components in the drawings are represented by the same reference numerals. It should be understood by those skilled in the art that the content specifically described below is illustrative rather than restrictive, and the protection scope of the present invention should not be limited thereto.

[0066] Based on the above technical problems, the inventor proposed a method, device and medium for stitching multispectral images of UAV remote sensing.

[0067] like Figure 1 As shown, one embodiment of the present invention proposes a method for stitching multispectral images of UAV remote sensing, including:

[0068] S10: Acquire N frames of single-band images and perform band registration correction to obtain N frames of multispectral images, wherein each frame of single-band image contains source images of M bands, wherein:

[0069] N and M are positive integers;

[0070] S20: performing band fusion on the N frames of multispectral images to obtain N frames of grayscale images, and performing splicing based on the N frames of grayscale images to obtain N local homography matrix sets;

[0071] S30: Through the nth local homography matrix set D n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band splicing results, and the M single-band splicing results are fused to obtain a preliminary splicing image, where:

[0072] n=1, 2,...,N; m=1, 2,...,M;

[0073] S40: Checking and correcting the line segment structure in the preliminary stitched image to obtain a final stitching result.

[0074] This scheme first performs band registration correction on the grayscale images of each band to ensure the consistency between the bands, and then obtains a set of local homography matrices by registering the fused single-band images to improve the efficiency of the algorithm. Finally, in view of the unique structural characteristics of farmland and road segments in agricultural scenes, line segment detection and correction are performed on the spliced ​​images to reduce the distortion and deformation that may be caused by splicing, and improve the accuracy and precision of the splicing results.

[0075] The multispectral camera of a drone is composed of several single-band cameras. Each camera has an independent optical path and records the data of the corresponding imaging band. Since the optical paths between the lenses are independent of each other, there are slight differences between the grayscale images of each band. Although these differences are relatively small, they will also affect the subsequent stitching. Therefore, it is necessary to perform band registration for the grayscale images of each band taken in the same scene.

[0076] like Figure 2 As shown in Figure 2, taking a UAV with 6 single-band spectral imagers as an example, first, the UAV collects N frames of video according to the preset route, such as Figure 3 As shown in the figure, the images of each band before band registration. For the single-band grayscale image data source of each band, there will be some differences between the small images of different bands taken in the same scene, so these images need to be corrected by band registration, such as Figure 4 As shown in the figure, the images of each band after band registration are shown. The corrected image ensures the spatial consistency between the bands, which is convenient for subsequent image stitching. Considering that the SIFT algorithm has scale invariance and can quickly and accurately complete grayscale image registration, the SIFT algorithm is used for band registration correction.

[0077] Specifically,

[0078] SIFT feature extraction is performed to obtain the feature descriptor of the image. The extracted feature descriptors are matched using the Euclidean distance algorithm. The transformation matrix is ​​calculated using the matching points. The image to be registered is transformed and interpolated to obtain the image after band registration.

[0079] In a possible implementation, the S10 includes:

[0080] S101: acquiring a video along a specific route through a multispectral camera of an unmanned aerial vehicle, wherein the video includes N frames of single-band images, and each frame of the single-band image includes source images of M bands;

[0081] S102: Using SIFT algorithm to perform band registration correction on the N frames of single-band images respectively:

[0082] Convolve the current single-band image with Gaussian kernel functions of different sizes, reduce the image after convolution at different levels to obtain the scale space, and perform extreme value detection on the source matrix and the network matrix in the scale space to obtain feature points, remove the edge response points in the feature point set, and assign direction values ​​to the key points;

[0083] Generate a feature point descriptor according to the direction value, and match the feature points to obtain matching points;

[0084] The matching points are used to calculate the transformation matrix to perform image transformation and difference on the source images of the M bands of the current single-band image, obtain the multispectral image after band registration, and update the source images of the M bands.

[0085] It should be understood that when using the SIFT algorithm for feature extraction, in order to ensure scale constancy, it is necessary to first create a scale (DOG) space and extract points of interest at different levels. In order to make the scale constant, the SIFT algorithm convolves the image with Gaussian kernel functions of different sizes, and then reduces the image after convolution at different levels to generate a DOG space.

[0086] Among them, the expression of the two-dimensional Gaussian kernel function G(x, y, σ) is shown in formula (1-1):

[0087]

[0088] Among them, σ is the standard deviation of the two-dimensional Gaussian kernel function, also called the scale factor. I(x,y) is the original image, and the scale space transformation formula is shown in formula (1-2):

[0089] L(x,y,σ)=I(x,y)*G(x,y,σ) (1-2)

[0090] The SIFT algorithm creates a Gaussian pyramid to obtain a thumbnail of an image. Based on downsampling, different models are used to apply Gaussian smoothing to each level of the image to obtain a Gaussian pyramid. In a Gaussian pyramid, the number of image groups (Octave) is the same as the number of layers (Interval), and multiple layers of images constitute a group. The first layer of the current image group is obtained by downsampling the third to last image of the previous group twice.

[0091] The SIFT algorithm uses the method of constructing DOG scale space to locate robust feature points. The DOG scale space is obtained by subtracting two adjacent images in the same layer. Then the function expression D(x, y, σ) for calculating the scale space is shown in formula (1-3):

[0092] D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y) (1-3)

[0093] Then, detection is performed in the DOG scale space to find feature points: the specific method to determine the position of feature points in the DOG scale space is to compare each point with all the adjacent points around it, including points on the same scale plane and points on the upper and lower planes, on each layer of each group of images. If the value of the pixel point is larger or smaller than the rest of the points, then the pixel point is considered to be a feature point. Many feature points are searched and found in the constructed scale space, but some feature points with edge effects need to be screened out. Because the curvature of the feature point is larger in the gradient direction at the boundary, and smaller in the orthogonal direction. The curvature of the selected feature point is proportional to the eigenvalue of the 2×2 Hessian matrix, and the H matrix is ​​shown in formula (1-4):

[0094]

[0095] Where D xx and D yy Denote the second-order partial derivatives along the x and y directions, respectively, D xy represents the second-order partial derivative along the x direction and then along the y direction. Let α and β be the eigenvalues ​​of the matrix H, and α>β, then the trace and eigenvalue of the matrix H are shown in equations (1-5) and (1-6):

[0096] Tr(H)=D xx +D yy =α+β (1-5)

[0097] Det(H)=D xx D yy -(D xy ) 2 =αβ (1-6)

[0098] The change of principal curvature is expressed by calculating the ratio of the trace and the eigenvalue. If you want to remove edge response points, you need to determine whether the difference between the point in the horizontal and vertical edge directions is large. If the difference is small, it means that this is not a point that needs to be removed, otherwise it needs to be removed. Therefore, r cannot exceed a certain threshold. As long as r is controlled to be less than the threshold, the ratio of the trace and the eigenvalue will also be less than the specified threshold, so that the difference in principal curvature in the two directions is within a certain range, and the edge response point is deleted.

[0099] Finally, generate the feature point descriptor: the feature descriptor contains information about the scale, position and direction of the feature point. Therefore, when generating the feature point descriptor, the first thing to calculate is the size and deviation angle of the gradient of the key point and the points around it; then use this information to draw the gradient statistical histogram, and the highest value in the histogram corresponds to the deviation angle of the key feature point. When calculating the direction of the feature point, the area around the feature point is divided into a 4×4 grid. Then, in these 16 small grids, each small grid will contain information about 8 directions, so there are a total of 128-dimensional vectors, which means that the feature descriptor has 128 dimensions. Then use the correct matching points to calculate the transformation matrix, perform image transformation and interpolation on the image to be registered, and obtain the image after band registration.

[0100] In a specific embodiment, the feature descriptor includes the scale, position and direction of the feature point.

[0101] In a specific embodiment,

[0102] The S20 includes:

[0103] S201: performing band fusion on the N multispectral images to obtain N grayscale images;

[0104] S202: performing single-band registration on the N grayscale images respectively;

[0105] S203: Adopting an approximate projection transformation algorithm to stitch the N grayscale images to obtain corresponding N local transformation matrix sets.

[0106] In a specific embodiment, the band-registered image is subjected to band fusion, and the fused grayscale image is subjected to grayscale image registration to obtain a local homography matrix set applicable to the entire band, which is then applied to the images of each band for image registration and fusion to obtain a spliced ​​image.

[0107] In the specific implementation process, the SIFT algorithm feature extraction and matching are first used to find the correct matching point pair, and then the linear transformation weights of the moving line are calculated in blocks, the source image is evenly divided into small grids, and the local projection transformation matrix is ​​calculated for each small grid in the image in a weighted manner for mapping to align with the target image. After obtaining the stitching result, the line segment structures such as farmland and roads in the image are detected and corrected to obtain a more accurate multispectral image stitching result.

[0108] In one possible implementation,

[0109] The S203 includes:

[0110] S2031: extracting and matching features of the N grayscale images respectively to obtain a pre-transformation matrix;

[0111] S2032: Perform block calculation of the linear transformation weight of the moving line, divide the grayscale image into blocks and evenly divide the image into small grids. The divided image is as follows: Figure 5 As shown, a weighted method is used to calculate the local projection transformation matrix for each small grid in the image for mapping; the moving line linear transformation algorithm obtains the local homography transformation matrix set D by calculating the weighted square sum error of the distance between the matching point in the current grid and the center of the grid. n ,in,

[0112] K is the total number of grids in the current grayscale image;

[0113] S2033: transforming the grayscale image using the local homography transformation matrix set, fusing and splicing the transformed image with the reference image using a weighted average fusion method,

[0114] in,

[0115] The local homography matrix of the kth grid of the nth grayscale image The calculation formula is shown in formula (1). The local transformation matrix is ​​estimated by solving the eigenvector corresponding to the minimum eigenvalue Values:

[0116]

[0117] Among them, a i is the linear mapping function between the matching points in the kth grid of the nth grayscale image, l is the number of matching points; d is the distance between the matching point and the center of the grid, is the weight corresponding to the feature point, W k is the diagonal matrix corresponding to the weight, and A is the mapping matrix form corresponding to the mapping function, where:

[0118]

[0119]

[0120] Among them, x k and k is the coordinate of the center point of the kth grid, and is the coordinate of the i-th matching point in the k-th grid, σ is the Gaussian scaling factor, which ranges between (0,1).

[0121] In a specific embodiment, Figure 2As shown, through the nth local homography matrix set D n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band spliced ​​images, and the M single-band spliced ​​images are fused to obtain a preliminary spliced ​​image. Specifically,

[0122] Through the first local homography matrix set D corresponding to the single-band image 1 1 Band registration is performed on the grayscale images of bands 1 to 6 of the first single-band image, and the second local homography matrix set D corresponding to the single-band image 2 is obtained. 2 Perform band registration on the grayscale images of bands 1 to 6 of the second single-band image; and so on;

[0123] The registered grayscale images of the same band are stitched together to obtain 6 single-band stitching results, which are then fused to obtain a preliminary stitching image.

[0124] In a specific embodiment, the S40 includes:

[0125] S401: Detect and extract line segment structures, extract line support regions in the preliminary spliced ​​image, and obtain line segment features through the line support regions, wherein:

[0126] The line support region is composed of pixel regions where line segments exist;

[0127] S402: Perform line segment correction: After obtaining the line segment features, use the least squares method to fit a line segment in the rectangular area, use the intersection of the line segments as feature points, sort their positional relationships with the center point of the rectangular area, calculate the perspective transformation matrix based on the correspondence between the transformed feature points and the original feature points, map all points in the area to be corrected, and complete the line segment correction.

[0128] In a specific embodiment, the S401 includes:

[0129] S4011: Use the four pixels to the lower right of each pixel to calculate the gradient corresponding to the pixel to obtain the position of the edge of the line segment. The gradient h of the pixel along the x direction and the y direction x (x,y) and h y The calculation formulas for (x, y) are shown in equations (4) and (5):

[0130]

[0131]

[0132] Where h(x, y) represents the grayscale value corresponding to the pixel with coordinates (x, y);

[0133] The calculation formula of gradient size H(x,y) and direction θ(x,y) is shown in formula (6):

[0134]

[0135] Among them, h x (x,y) and h y (x, y) are the gradients of the pixel along the x and y directions respectively;

[0136] S4012: Perform pseudo gradient sorting: divide the gradient value into 1024 equal parts from large to small, and take out a seed pixel from each segment for sorting;

[0137] Calculate the direction of the candidate rectangular area formed by the pixel points, and determine whether the pixel point belongs to the rectangular area by comparing the direction of the candidate pixel point and the candidate rectangular area, thereby completing the rectangular area division;

[0138] S4013: performing approximate calculation of rectangular areas: calculating the approximate rectangular areas corresponding to the line segments, the size of the rectangular areas covering the entire line support area, and each rectangular area represents a candidate line segment;

[0139] The center point c of the rectangular area x The calculation formula is shown in formula (7):

[0140]

[0141] Where (c x ,c y ) represents the coordinates of the center point, Region represents the current region, G(j) represents the gradient value of pixel j in the region, subscript j is used to traverse all pixels in the region, x(j) and y(j) represent the horizontal and vertical coordinates of pixel j respectively;

[0142] S4014: Improve the rectangular area and select the correct line support area: calculate the number of false alarms NFA value in the rectangular area. If the calculation result is less than a preset threshold, the preset threshold is, for example, 1, then the current rectangular area is considered to be the correct line support area. The calculation formula of the number of false alarms NFA value in the rectangular area r is shown in formula (8):

[0143] NFA(r,t)=(PQ) 2.5 γB(f(r),k(r,t),p) (8)

[0144] Among them, P and Q are the image resolution sizes, B() represents the binomial distribution, f(r) is the number of pixels in the rectangular area r, k(r,t) represents the number of points within the class in the rectangular area r, p is the probability of gradient angle alignment, and γ is the possible number of different values ​​of p.

[0145] It should be noted that gradient angle alignment means that all horizontal line direction angles in the rectangular area are considered aligned if the deviation between the rectangular direction angles is less than a threshold, and the threshold is, for example, 5°, 6° or 4°.

[0146] like Figure 6 As shown, a schematic diagram of the structure of a computer device provided by the second embodiment of the present invention. A computer device suitable for implementing the UAV remote sensing multispectral image stitching method provided in the above embodiment includes a central processing module (CPU), which can perform various appropriate actions and processes according to a program stored in a read-only memory (ROM) or a program loaded from a storage part to a random access memory (RAM). Various programs and data required for the operation of the computer device are also stored in the RAM. The CPU, ROM and RAM are connected to each other via a bus. The input / output (I / O) interface is also connected to the bus.

[0147] The following components are connected to the I / O interface: an input part including a keyboard, a mouse, etc.; an output part including a liquid crystal display (LCD), etc. and a speaker, etc.; a storage part including a hard disk, etc.; and a communication part including a network interface card such as a LAN card, a modem, etc. The communication part performs communication processing via a network such as the Internet. A drive is also connected to the I / O interface as needed. Removable media, such as a magnetic disk, an optical disk, a magneto-optical disk, a semiconductor memory, etc., are installed on the drive as needed so that the computer program read therefrom is installed into the storage part as needed.

[0148] In particular, according to this embodiment, the process described in the flowchart above can be implemented as a computer software program. For example, this embodiment includes a computer program product, which includes a computer program tangibly contained on a computer-readable medium, and the computer program includes program code for executing the method shown in the flowchart. In such an embodiment, the computer program can be downloaded and installed from a network through a communication part, and / or installed from a removable medium.

[0149] The flow chart and schematic diagram in the accompanying drawings illustrate the possible architecture, function and operation of the system, method and computer program product of the present embodiment. In this regard, each square box in the flow chart or schematic diagram can represent a part of a module, program segment or code, and a part of the above-mentioned module, program segment or code contains one or more executable instructions for realizing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the square box can also occur in a sequence different from that marked in the accompanying drawings. For example, two square boxes represented in succession can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each square box in the schematic diagram and / or flow chart, and the combination of the square boxes in the schematic diagram and / or flow chart can be implemented with a dedicated hardware-based system that performs a specified function or operation, or can be implemented with a combination of dedicated hardware and computer instructions.

[0150] A fourth embodiment of the present application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements:

[0151] S10: Acquire N frames of single-band images and perform band registration correction to obtain N frames of multispectral images, wherein each frame of single-band image contains source images of M bands, wherein:

[0152] N and M are positive integers;

[0153] S20: performing band fusion on the N multispectral images to obtain N grayscale images, and performing splicing based on the N grayscale images to obtain N local homography matrix sets;

[0154] S30: Through the nth local homography matrix set D n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band spliced ​​images, and the M single-band spliced ​​images are fused to obtain a preliminary spliced ​​image, wherein,

[0155] n=1, 2,...,N; m=1, 2,...,M;

[0156] S40: Checking and correcting the line segment structure in the preliminary stitched image to obtain a final stitching result.

[0157] In practical applications, the computer-readable storage medium can adopt any combination of one or more computer-readable media. The computer-readable medium can be a computer-readable signal medium or a computer-readable storage medium. The computer-readable storage medium can be, for example, but not limited to, a system, device or device of electricity, magnetism, light, electromagnetic, infrared, or semiconductor, or any combination of the above. More specific examples (non-exhaustive list) of computer-readable storage media include: an electrical connection with one or more wires, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above. In the present embodiment, the computer-readable storage medium can be any tangible medium containing or storing a program, which can be used by an instruction execution system, device or device or used in combination with it.

[0158] Computer-readable signal media may include data signals propagated in baseband or as part of a carrier wave, which carry computer-readable program code. Such propagated data signals may take a variety of forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination of the above. Computer-readable signal media may also be any computer-readable medium other than a computer-readable storage medium, which may send, propagate, or transmit a program for use by or in conjunction with an instruction execution system, apparatus, or device.

[0159] The program code embodied on the computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.

[0160] Computer program code for performing the operation of the present application can be written in one or more programming languages ​​or a combination thereof, including object-oriented programming languages, such as Java, Smalltalk, C++, and conventional procedural programming languages, such as "C" language or similar programming languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a separate software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In the case of a remote computer, the remote computer can be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or can be connected to an external computer (e.g., using an Internet service provider to connect through the Internet).

[0161] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not limitations on the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the embodiments here. All obvious changes or modifications derived from the technical solution of the present invention are still within the protection scope of the present invention.

Claims

1. A method for stitching multispectral images of UAV remote sensing, characterized in that: include: S10: Acquire N frames of single-band images and perform band registration correction to obtain N frames of multispectral images, wherein each frame of single-band image contains source images of M bands, wherein: N and M are positive integers; S20: performing band fusion on the N frames of multispectral images to obtain N grayscale images, and performing splicing based on the N grayscale images to obtain N local homography matrix sets; S30: Through the nth local homography matrix set D n The grayscale images of the mth band of the nth single-band image are respectively registered, and the registered images of the same band are spliced ​​to obtain M single-band splicing results, and the M single-band splicing results are fused to obtain a preliminary splicing image, where: n=1, 2,...,N; m=1, 2,...,M; S40: Checking and correcting the line segment structure in the preliminary stitched image to obtain a final stitching result.

2. The method according to claim 1, characterized in that The S10 includes: S101: acquiring a video along a specific route through a multispectral camera of an unmanned aerial vehicle, wherein the video includes N frames of single-band images, and each frame of the single-band image includes source images of M bands; S102: performing band registration correction on the N frames of single-band images respectively by using the SIFT algorithm: performing convolution with the current single-band image by using Gaussian kernel functions of different sizes, reducing the image after the convolution at different levels to obtain a scale space, and performing extreme value detection on the source matrix and the mesh matrix in the scale space to obtain feature points, while removing edge response points in the feature point concentration, and assigning direction values ​​to key points; Generate a feature point descriptor according to the direction value, and match the feature points to obtain matching points; The matching points are used to calculate the transformation matrix to perform image transformation and difference on the source images of the M bands of the current single-band image, obtain the multispectral image after band registration, and update the source images of the M bands.

3. The method according to claim 2, characterized in that The feature point descriptor includes the scale, position and direction of the feature point.

4. The method according to claim 1, characterized in that: The S20 includes: S201: performing band fusion on the N frames of multispectral images to obtain N grayscale images; S202: performing single-band registration on the N grayscale images respectively; S203: Adopting an approximate projection transformation algorithm to stitch the N grayscale images to obtain corresponding N local transformation matrix sets.

5. The method according to claim 4, characterized in that The S203 includes: S2031: extracting and matching features of the N grayscale images respectively to obtain a pre-transformation matrix; S2032: Divide the grayscale image into small grids evenly. The moving line linear transformation algorithm calculates the weighted square sum of the distances between the matching points in the current grid and the center of the grid to obtain the local homography transformation matrix set D. n ,in, K is the total number of grids in the current grayscale image; It refers to the local homography matrix of the kth grid of the nth grayscale image; S2033: transforming the grayscale image using the local homography transformation matrix set, and fusing and splicing the transformed image with the reference image using a weighted average fusion method.

6. The method according to claim 5, characterized in that The local homography matrix of the kth grid of the nth grayscale image The calculation formula is shown in formula (1). The local transformation matrix is ​​estimated by solving the eigenvector corresponding to the minimum eigenvalue Values: Among them, a i is the linear mapping function between matching points, l is the number of matching points; d is the distance between the matching point and the center of the grid, is the weight corresponding to the feature point, W k is the diagonal matrix corresponding to the weight, and A is the mapping matrix form corresponding to the mapping function, where: Among them, x k and k is the coordinate of the center point of the kth grid, and is the coordinate of the i-th matching point in the k-th grid, σ is the Gaussian scaling factor, which ranges between (0,1).

7. The method according to claim 1, characterized in that The S40 includes: S401: Detect and extract line segment structures, extract line support regions in the preliminary spliced ​​image, and obtain line segment features through the line support regions, wherein: The line support region is composed of pixel regions where line segments exist; S402: Perform line segment correction: Use the least squares method to fit a line segment in the rectangular area, use the intersection of the line segments as feature points, sort their positional relationships with the center point of the rectangular area, calculate the perspective transformation matrix based on the correspondence between the transformed feature points and the original feature points, map all points in the area to be corrected, and complete the line segment correction.

8. The method according to claim 7, characterized in that The S401 includes: S4011: Use the four pixels to the lower right of each pixel to calculate the gradient corresponding to the pixel to obtain the position of the edge of the line segment. The gradient h of the pixel along the x direction and the y direction x (x,y) and h y The calculation formulas for (x, y) are shown in equations (4) and (5): Where h(x, y) represents the grayscale value corresponding to the pixel with coordinates (x, y); The calculation formula of gradient size H(x,y) and direction θ(x,y) is shown in formula (6): Among them, h x (x,y) and h y (x, y) are the gradients of the pixel along the x and y directions respectively; S4012: Perform pseudo gradient sorting: divide the gradient value into 1024 equal parts from large to small, and take out a seed pixel from each segment for sorting; Calculate the direction of the candidate rectangular area formed by the pixel points, and determine whether the pixel point belongs to the rectangular area by comparing the direction of the candidate pixel point and the candidate rectangular area, thereby completing the rectangular area division; S4013: performing approximate calculation of rectangular areas: calculating the approximate rectangular areas corresponding to the line segments, the size of the rectangular areas covering the entire line support area, and each rectangular area represents a candidate line segment; The center point c of the rectangular area x The calculation formula is shown in formula (7): Where (c x ,c y ) represents the coordinates of the center point, Region represents the current region, G(j) represents the gradient value of pixel j in the region, subscript j is used to traverse all pixels in the region, x(j) and y(j) represent the horizontal and vertical coordinates of pixel j respectively; S4014: Improve the rectangular area and select the correct line support area: Calculate the number of false positives NFA value in the rectangular area. If the calculated result is less than a preset threshold, the current rectangular area is considered to be the correct line support area. The calculation formula of the NFA value is shown in formula (8): NFA(r,t)=(PQ) 2.5 γB(f(r),k(r,t),p) (8) Among them, P and Q are the image resolution sizes, B() represents the binomial distribution, f(r) is the number of pixels in the rectangular area r, k(r,t) represents the number of points within the class in the rectangular area r, p is the probability of gradient angle alignment, and γ is the possible number of different values ​​of p.

9. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the method according to any one of claims 1 to 8 is implemented.

10. A computer storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the method according to any one of claims 1 to 8 is implemented.

Images