AUTOMATED REGISTRATION OF FULL MOON REMOTE SENSING IMAGES BASED ON TRIANGULATED NETWORK CONSTRAINTS

: The registration of full-moon remote sensing images constitutes a pivotal stage in the fusion analysis of multiple lunar remote sensing datasets. Addressing prevailing issues in automatic registration, such as the broad width of full-moon data, signiﬁcant internal distortion, and texture distortion in high-latitude regions, this paper proposes a method for automatic matching and correction based on triangulation constraints. The approach employs a matching strategy progressing from coarse to ﬁne and from sparse to dense. It optimizes and combines multiple existing matching algorithms, enhances the extraction of initial network points, constructs irregular triangulation networks using these points, conducts dense matching with each triangulation network as a basic unit, and introduces a geometric correction method based on triangulation network + grid (TIN + GRID) for the registration of full-moon data. For the matching of full-moon remote sensing images in high-latitude regions, a novel approach involving memory projection forward transformation-matching-projection inverse transformation is adopted. Through registration experiments with full-moon image data and an analysis of registration accuracy at different latitudes, the average mean square error is found to be less than 2 pixels. These results signify the efﬁcacy of the proposed method in effectively addressing the automatic registration challenges encountered in full-moon remote sensing images.


INTRODUCTION
Full-moon image data serves as the foundational dataset for scientific investigations into lunar topography, mineral exploration, and geological evolution.Since the inception of lunar exploration missions in the 1950s, a substantial volume of lunar image data has been acquired [Tong, 2022, Xu, 2022].The Clementine mission produced a 100-meter resolution full-moon stitched map [Scholten, 2012, Robinson, 2010b]; Barker et al. produced a DOM product with a resolution of 59 meters covering latitudes of 65°north and south based on the images of the Japanese "Moon Goddess" satellite [Barker, 2016a]; the Indian "Chandrayaan-1" satellite obtained a large amount of images with a resolution of 5 meters [Goswami, 2009]; the Lunar Reconnaissance Orbiter (LRO) launched by NASA obtained images covering almost the entire moon [Robinson, 2010a, Jha, 2019].CE-1 and CE-2 also obtained full-moon orthophotomap images with different resolutions such as 120 meters and 7 meters [Li, 2013, Li, 2018].Although different lunar probes have obtained a large amount of remote sensing data on the lunar surface, the control data and reference benchmarks used in the mapping products produced by different missions are different, resulting in geometric deviations between lunar DOM products [Di, 2019,Wu, 2011,Wu, 2014].This severely restricts the comprehensive application of multiple mission data for scientific research.Therefore, it is necessary to carry out the registration research of different types of lunar DOM products and unify various data into the same reference benchmark.
Many scholars have conducted research on lunar image registration methods [Xin, 2018, Liu, 2023, Yang, 2020, Yang, 2023, Bo, 2022, Jiang, 2021], but there are few studies on the registration of full moon images.Barker et al. used a block registration method to produce a full moon DEM; in this method, blocks are independent from each other and are not truly globally corrected [Barker, 2016b]; Bo et al. proposed an automatic registration method for lunar data based on spherical triangular grids; in this method, feature point extraction and matching are cut out for the north and south poles for separate processing, which is not truly based on a full moon image registration.The registration of full moon images faces many difficulties.First, full moon data are spliced from a large number of single-scene images taken at different times, and the geometric accuracy of different regions within the data is inconsistent.Second, in the spherical coordinate system, there is a large amount of redundant storage in images of high-latitude regions, and there is severe texture distortion visually, making it difficult to accurately extract the corresponding points required for registration.Third, how to ensure the accuracy and internal accuracy consistency of the corrected full moon image is also a difficulty faced by registration.
Common image matching methods mainly include region-based methods, feature-based methods, and deep learning-based methods.The region-based automatic matching method is a method that matches images based on a certain template size, using the center of the template window as the point of correspondence, and implementing image matching based on a certain feature similarity measure.Common methods include correlation coefficient and mutual information [Suri, 2010, Gutiérrez-Becker, 2016].For images with inconsistent internal accuracy, a single size template may not match points.The feature-based method is to extract features between images and use the similarity between features for matching [Ye, 2014, Li, 2016, Li, 2015, Yang, 2021]; common feature-based methods include methods based on point, line, and surface features [Yu, 2008, Sui, 2015, Zhang, 2007].The full-moon image is composed of many images from different periods, and there are nonlinear radiometric differences within it.Conventional feature matching methods may have a high mismatch rate, and even the number of correctly matched points is far less than the number of mismatched points.The deep learning-based matching method is to automatically select image features and determine the correspondence between features through learning [Hao, 2023, Ma, 2021, Bürgmann, 2019].This method requires extensive training of matching samples to be applied, and currently lacks samples for lunar image registration.
To address these challenges, this paper proposes an automatic registration method for full-moon remote sensing images based on triangular network constraints.The method employs a matching strategy progressing from coarse to fine and from sparse to dense, optimizing and amalgamating multiple existing matching algorithms, refining the extraction of initial network points, and constructing irregular triangular networks using these points.Subsequently, dense matching is conducted with each triangular network as a fundamental unit.Finally, geometric correction is implemented using the triangular network + grid method.Throughout the matching process, techniques such as dynamic memory projection transformation in high-latitude regions are utilized to resolve registration issues associated with internal accuracy inconsistencies and texture distortions in such regions for full-moon data.The primary contributions of this paper are outlined as follows: 1)To address texture distortion and matching challenges in highlatitude regions of full-moon remote sensing images, the method adopts the concept of memory projection forward transformationmatching-projection inverse transformation.This ensures that the entire matching process remains unaffected by latitude-based texture variations.
2)The paper introduces an automatic image matching and correction method based on triangulation constraints.Employing a coarse-to-fine, sparse-to-dense matching strategy and optimizing a combination of existing matching algorithms, it resolves the registration challenges associated with full-moon data, ensuring high registration accuracy.
3)The proposed method includes a technique for constructing a corrected grid (GRID) on top of the Triangulated Irregular Network (TIN) established at control points.This enhances orthoprojection output efficiency, mitigates the impact of a single erroneous control point causing significant local misalignment, and addresses internal consistency issues in the corrected fullmoon image.

Technical route
In order to solve the problem of high-precision registration of data throughout the month, this article proposes a registration method based on triangular network constraints.The specific idea is as follows: first, match a certain number of corresponding point on the images to be matched and the reference image, which are called the initial matching points; then construct irregular Delaunay triangular networks based on the initial matching points; then perform dense matching under the constraints of each triangular network; finally, based on the final obtained corresponding points, use the method of triangular network + grid (TIN + GRID) for geometric correction.The registration process is shown in Figure 1 [Harris, 1988].This algorithm mainly detects corners by calculating the curvature and gradient of the point.
Harris feature point extraction algorithm is not only computationally simple, but also has high stability and robustness, and can accurately detect feature points under image rotation, grayscale changes, and noise interference.
In order to extract evenly distributed feature points on the fullmonth data, the Harris operator is used.The Harris operator extracts feature points by setting a specific threshold.If a threshold is set on the entire image to extract feature points, it will result in uneven distribution of feature points.The block Harris operator divides the image into regular image blocks by latitude and longitude; each block uses the Harris operator to calculate the feature value, and selects the point with the largest feature value in the current block as the feature point;this avoids the uneven distribution of feature points, as shown in Figure 2 below.According to Equation 1, It can be calculated based on the aspect ratio, so that the starting position and step size of each image block can be determined.The image pyramid describes and records the grayscale information of images at different levels.High-resolution images are placed at the lower levels, and low-resolution images are placed at the upper levels, forming a pyramid shape.This article uses the coarse-to-fine features of the pyramid to accelerate the efficiency of matching.First, automatic matching is performed at the top level of the pyramid, and the matched points are used to optimize the initial matching position of the next level of the pyramid.Matching is performed level by level until the original layer of image data is reached.
If the two-scene full-month remote sensing images are directly matched point by point with a large search radius R, it will be very time-consuming.Therefore, in this paper, a certain number of corresponding points are automatically matched in the pre-and post-processing remote sensing images, which we call the initial matching points.The matching of the initial matching points can use a larger search radius to deal with situations where the initial accuracy deviation is large.The accuracy of the initial matching points has a crucial impact on subsequent processing, so it should be improved as much as possible to ensure that the initial matching points are correct.
There are many methods for automatic matching of image corresponding points, including the commonly used correlation coefficient matching method [Zitová, 2003] and the channel features of orientated (CFOG) that supports automatic matching of heterogeneous images [Ye, 2019].The correlation coefficient matching method is suitable for automatic matching of homologous remote sensing images and has high efficiency, while the CFOG can be used for automatic matching of homologous and heterogeneous remote sensing images, with high matching success rate and low efficiency and many error points.
The generation of initial matching points in this article adopts a matching method based on correlation coefficient and CFOG.Each feature point must satisfy both matching methods simultaneously to be finally identified as a successful match, thus improving the accuracy of initial matching points.The search radius R of matching parameters can be set to 200, and the feature points are evenly distributed in the overlapping area of the preand post-remote sensing images.The number of feature points N can be determined based on the image size.The more feature points, the slower the efficiency, and the less feature points, the less dense the constructed triangular network.Through experiments, it is found that usually 10,000 feature points can be extracted per 100 million pixels to set up an equal proportion.
The correlation coefficient and CFOG matching algorithm are described as follows: The correlation coefficient matching method uses the correlation coefficient (standardized covariance) as a similarity measure.In statistics, the correlation coefficient is used to represent the correlation between two random variables, and in image matching, it can be used to represent the similarity between two images of the same size.In practice, it is not necessary to compute the directional gradient channel go separately for each layer; the directional gradient channel eigenvalues for each layer can be computed by using the gradient magnitudes (gx, gy) in the horizontal and vertical directions: where θ = the gradient direction of the division the purpose of the absolute value is to limit the gradient direction to [0, π], which can better deal with the case of gradient reversal among multimodal remote sensing images.
(2) y constructing a three-dimensional Gaussian convolution kernel, the three-dimensional multi-directional gradient feature maps constructed in the previous step are convolved to finally form CFOG descriptors.Among them, the 3D Gaussian convolution kernel is not a strictly 3D Gaussian function, but a 2D Gaussian kernel in the X and Y directions and a kernel in the gradient direction.
where * = the the convolution σ = the standard deviation of the Gaussian convolution kernel σ is the standard deviation of the Gaussian convolution kernel, which is not strictly a three-dimensional Gaussian function in three-dimensional space, but a two-dimensional Gaussian kernel in the X and Y directions and a kernel [121] T in the gradient direction (hereinafter referred to as the Z direction).
The gradient in the Z-direction is smoothed by convolution in the Z-direction, which reduces the effect of directional distortion between the two images due to local geometrical deformations and changes in gray intensity; the final directional gradient channel feature is formed by normalization operation on g o o .
2) Frequency domain spatial matching: Fourier transform is used to transform the feature template from the spatial domain to the frequency domain, and phase correlation is used as the similarity measure for template matching.

Constructing a spherical Delaunay triangulation
An irregular triangular mesh is a series of discrete data points connected into a series of continuous triangular meshes, the size and shape of which depend on the location and density of the discrete data points.Using the discrete data obtained from all sampling points, according to the principle of optimal combination, these discrete points (the vertices of the triangles) are connected into a continuous triangular surface (in the connection, as far as possible, to ensure that each triangle is an acute triangle or the length of the three sides of the triangle is approximately equal), the shape and size of the triangular surface depends on the location and density of the irregularly distributed measurement points or nodes.
In order to be able to evenly distribute the encrypted points on top of the initial distribution points, Delaunay triangulation was used to construct the full moon triangulation control network.Due to the problems such as the difficulty of processing the boundary of planar Delaunay triangulation network and the unevenness of triangulation network on the boundary [Zhou, 2007], spherical Delaunay triangulation network is used to construct the network in this paper.
The main algorithms applicable to the spherical surface to generate Delaunay triangular mesh are Recursive Subdivision Algorithm, Spherical Centroid Method, Insertion Algorithm, Centroid Projection Method and so on.Centroid Projection Method) and so on [Fang, 1995, Shi, 2006, Yang, 2018]; among them, the mesh generated by the Centroid Projection Method has a better shape uniformity, and this method is chosen for the spherical Delaunay triangular mesh construction in this paper.

Triangulation constrained dense matching
Based on the irregular Delaunay triangulation constructed at the initial matching point, dense matching is performed under the constraints of each triangulation.The matching method still uses the pyramid-level image matching strategy, and the matching operator uses the CFOG matching algorithm.Each triangular network consists of three initial matching points, and the three points can be used to solve the parameters of the affine transformation.The affine transformation model formula is as follows: The above formula describes the affine transformation relationship between two-dimensional points x and y and f (x) and f (y).The six parameters a0, a1, a2, b0, b1, b2 of the affine transformation are the objects of the equation system solution.
Due to the established coordinate affine transformation relationship between the images to be matched and the reference images at the initial matching point, there are several points to note when using the dense matching method within the triangulation network: (1)Due to the high accuracy of the three initial matching points in the triangulation network, it can provide more accurate location prediction.Search radius can be set very small to speed up matching efficiency; (2) The higher the degree of intensive matching, the longer the matching time; a comprehensive value can be calculated based on the matching requirements and matching time; (3) After the triangulation network is densely matched, the error values of each densely matched point should be calculated according to the affine parameters formed by the initial matching points, and the erroneous matching points should be deleted according to the error threshold.

Special treatment in high latitude regions
The fullmonth data used in this article is in the Moon2000 lunar geographical coordinate system, which is a spherical coordinate system.In the spherical coordinate system, images in highlatitude regions have a large amount of redundant storage, resulting in severe visual distortions that are not conducive to feature point extraction and matching.
To solve this problem, before matching the initial matching points, it is first determined whether the latitude of the matching point is higher than 60°.If it is higher than 60°, the image within a certain radius around the matching point is dynamically converted to the polar location projection, and the homonym matching is performed under the polar location projection.Considering the large amount of computation required for point-by-point dynamic projection, in the dense matching step constrained by the triangular network, it is first calculated whether there are any points with latitudes exceeding 60°among the three initial matching points that form the triangular net- This method proposes a method of constructing a corrected grid (GRID) based on the construction of a TIN using control points, aiming to improve the efficiency of the algorithm while reducing the impact of an erroneous control point that leads to significant local misalignment.
First, calculate the six affine transformation parameters for each triangle using formula ( 6  Point P is between pixel points 11, 12, 21, and 22, and the grayscale resampled value I(P) is where I11, I12, I21, I22 = the gray values of the four adjacent pixels of P x, y = the coordinates of P, usually in decimal ∆x and ∆y are as follows, which are the differences from their integer values.

Experiment data
The full-moon registration method designed in this article was verified using WAC data [WAC., 2013] and CE-1 data [CE1., 2020]; WAC data was used as the reference image.WAC data was downloaded from the USGS website, and CE-1 data was downloaded from the Lunar and Planetary Data Release System website.WAC data is a full-moon DOM data with a resolution of 100m/pixel, which was created by stitching more than 15,000 images obtained by the Wide Angle Camera

Analysis of matching results and matching accuracy
Matching parameters: refer to the method of million standard sub-frames for blocking, where the latitudes between 0 degrees and 60 degrees are blocked according to the longitude difference of 6 degrees and the latitude difference of 4 degrees, the latitudes between 60 degrees and 76 degrees are blocked according to the longitude difference of 12 degrees and the latitude difference of 4 degrees, the latitudes between 76 degrees and 80 degrees are blocked according to the longitude difference of 24 degrees and the latitude difference of 4 degrees, and the latitudes above 80 degrees are blocked as one block.The whole moon is divided into 2072 blocks; the number of seed points per block is 10; the initial search radius for matching is set to 30 pixels.
After using correlation coefficient and CFOG algorithm matching, as well as RANSAC to remove erroneous points, a total of 4,120 initial matching points were obtained.Figure 6 shows the distribution of matching results in some low-and mid-latitude regions and high-latitude regions.700 checkpoints were evenly distributed to check the matching accuracy.The maximum error was 1.2 pixels, the minimum error was 0.5 pixels, the average error was 0.96 pixels, and the median error was 0.99 pixels.This indicates that the initial matching points have achieved pixel-level accuracy, providing a reliable accuracy guarantee for subsequent dense matching.

X direction (pixels)
Y direction (pixels) XY direction (pixels) Based on the initial point, a total of 4118 Delaunay triangulation networks were generated using spherical generation.Figure 7 shows the distribution of some triangulation networks in lowand mid-latitude regions and high-latitude regions.
Based on the irregular Delaunay triangulation constructed at the initial matching point, dense matching is performed under the constraints of each triangulation.Matching parameters: using the block Harris algorithm, the initial search radius is 10 pixels, and 100,000 matching points are arranged overall.The triangulation of each matching point is calculated, and the error value of each dense matching point is calculated using the affine parameters calculated from the three initial matching points of the triangulation.The error threshold is used to delete erroneous matching points, resulting in a total of 68,213 encrypted points.Uniformly extract 100 triangulation network checkpoints to check the matching accuracy of the encrypted points.The maximum value is 2.2 pixels, the minimum value is 0.36 pixels, the average value is 1.58 pixels, and the mean square error is 1.60 pixels.This indicates that the accuracy of the encrypted points meets the matching requirements.

X direction (pixels)
Y direction (pixels) XY direction (pixels)  3.3 Analysis of the whole-month registration results

Qualitative analysis
This section mainly focuses on visual inspection, extracting 110 regions evenly at certain intervals from the monthly data, and manually inspecting the registration effect.Figure 8 shows the distribution of inspection regions and the inspection results for some low-latitude, midlatitude, and high-latitude regions.
As can be seen from the image rolling curtain diagram, there are deviations in different latitudes before registration, and the deviations are not consistent in different regions.However, after registration, the deviations are smaller, indicating that this method is suitable for large-scale, internally distorted image registration.From low-latitude to high-latitude regions, the superimposition between the registered image and the reference image is very good, reaching a precision of 1.5 pixels, indicating that this method solves the registration problem at high latitudes and improves the internal consistency of the corrected full-moon image.The error calculation formula for each checkpoint is: The final calculation formula for the mean square error is:  To further analyze the deviation of different latitudes before and after registration, a distribution diagram of errors with latitude was created.The specific method is to take ten degrees as an interval from 0-90 degrees, and calculate the average value of errors within each interval.
From the distribution map of different latitude errors, it can be seen that before registration, the deviation between different latitudes is quite large, with the largest deviation occurring in the polar regions.After registration using the method proposed in this paper, the deviation between different latitudes tends to be consistent, indicating that the internal distortion of the registered images is relatively small.As shown in Figure 10: In addition, in order to verify whether the method proposed in this article is applicable to high-latitude regions, the area above 75°latitude in the northern hemisphere before and after registration was cropped and converted to the polar stereographic projection.80 checkpoints were evenly distributed at certain intervals, and the error values before and after registration were calculated.From the comparison of errors, the geometric position accuracy of high-latitude regions has been greatly improved after registration, which reflects the effectiveness of the high-latitude processing method adopted in this article.As shown in Figure 11, Table 4:

CONCLUSION
In view of the challenges of fully automatic registration such as large texture differences and internal distortions in full-month data, this paper proposes an automatic image registration method based on triangular network constraints.The method first performs fine-grained automatic matching to obtain initial matching points.Using the initial matching points, an irregular triangular network is established.Taking each triangle as a unit, the location constraints of the three vertices are used to quickly perform dense matching on the feature points inside the triangle and calculate gross errors.For the initial matching points and encrypted points in high-latitude regions, a single matching point and a single triangular network are used as units, respectively.The idea of memory projection forward transformationmatching-projection inverse transformation is adopted to solve the problem of texture distortion matching in high-latitude regions.This matching method not only improves efficiency, but also ensures the distribution uniformity of matching points in full-month image data.By using the TIN+GRID geometric correction method, high-precision registration results are ensured, especially solving the problem of local error points affecting the overall registration result.Finally, by qualitatively and quantitatively analyzing the results of full-month registration, the effectiveness of this method is proven.

Figure 1 .
Figure 1.Overall technical scheme for registration Figure 2. Comparison of the effect of Harris operator blocking and non-blocking ) where R(x, y) = the correlation coefficient of the two pieces of image data E(X), E(Y ) = the mean gray value of the two pieces of image data D(X), D(Y ) = the variance of the two pieces of image data E(XY ) = the mean value of the corresponding points of the two pieces of image data after multiplication |R(x, y)| = The closer the correlation coefficientE(XY ) is the mean value of the corresponding points of the two pieces of image data after multiplication, and their definitions are the same as those defined in general statistical theory.The correlation coefficient indicates the degree of similarity of the linear relationship between images X and Y.The closer the correlation coefficient |R(x, y)| is to 1 or -1, the more obvious the degree of linear similarity between the images.When the correlation coefficient is greater than the set threshold, the feature point matching is considered successful.The full moon data produced by different sensors often have nonlinear radiometric differences, which can cause the correlation coefficient algorithm to fail in multi-source image feature matching.The CFOG (Chanel Feature of Orientated Gradient) matching algorithm is a pixel-by-pixel feature representation of the image by calculating the histogram of gradient direction for each pixel.The CFOG matching algorithm mainly consists of two processes: establishing CFOG descriptors and frequency domain spatial matching.1)Constructing the CFOG descriptor(1) Directional gradient channel calculation: Calculate the gradient channel values in multiple directions for a given image, and arrange all the gradient channel values in the Z direction to form a three-dimensional multi-directional gradient feature map.Define the directional gradient channel g o for each layer as: direction of division = the result when the value is positive, otherwise it is 0 Figure 3.Comparison of high-latitude geographical images and planar images ); then, define the size G of the corrected grid, which can be adjusted adaptively based on the output size of the image.Usually, G is set to 10. Divide the corrected image coordinate range into G*G squares.Assuming that the output image width is W and height is H, there are W/G * H/G squares.Determine which triangle the center point (x1, y1) of each corrected grid falls within, and then use the affine transformation parameters of the triangle to calculate the pixel coordinate value (x2, y2) of the grid center point in the image before correction.Finally, we obtain the corrected image coordinates of the G*G grid center points and the image coordinates before correction, which are stored in two memory twodimensional arrays Grid x[W/G, H/G] and Grid y[W/G, H/G].According to the corrected grid established in the previous step, the coordinate mapping relationship between the corrected image and the image to be corrected can be quickly located.Cycle each pixel P (x1, y1) of the corrected image, then obtain the grid coordinate P grid (x grid, y grid) of the floating point number by the coordinate of each pixel is divided the grid size G.Take the Grid x[W/G, H/G] and Grid y[W/G, H/G] coordinate values of the four nearest grid points to point P grid, and calculate the value obtained by bilinear interpolation algorithm, which should be taken from the coordinate value P src(x2, y2) of the rectified image before.The four pixel values around the P src point are bilinearly interpolated to obtain the final coordinate value of P (x1, y1).This step uses bilinear interpolation twice.The first bilinear interpolation is for coordinates, and the second bilinear interpolation is for grayscale values.

Figure 4 .
Figure 4. Schematic diagram of bilinear interpolation algorithm Figure 5. Full-month registration data

Figure 7 .
Figure 7.The effect of spherical Delaunay triangulation at different latitudes

Figure 8 .
Figure 8. Visual inspection area distribution and comparison of partial areas before and after registration (left is the reference image)

Figure 9 .
Figure 9. Distribution of quantitative inspection checkpoints checkpoints in longitude direction between the matching image and the reference image ∆y = the difference of checkpoints in the latitude direction between the matching image and the reference image ∆xy = the comprehensive error for a single checkpoint ∆XY = the comprehensive mean square error of single scene image n = the number of checkpointsThis study separately calculated the deviations between the registration image and the reference image before registration, and the deviations between the registered image and the reference image after registration.From the value of the mean square error before and after registration, it can be seen that the geometric position accuracy has been greatly improved after re-

Figure 10 .
Figure 10.Distribution of different latitude errors before and after registration

Figure 11 .
Figure 11.Distribution of checkpoints in the northern hemisphere (above 75 degrees)

Table 1 .
Matching accuracy at the initial matching point

Table 2 .
Matching accuracy of encryption pointst

Table 3 .
Error before and after registration

Table 4 .
Error before and after registration in high latitude areas