A method and system for image mosaic fusion based on topology analysis

Abstract

The application discloses a kind of based on topological analysis's graphic image splicing fusion method and system, belong to image processing field, the method includes: obtaining first image and second image, generate first binary image and second binary image;Respectively topological feature analysis is carried out, and the topological skeleton of image foreground area is extracted and topological invariants is calculated, and initial matching point pair set is generated;After alignment, the first image is generated;Image splicing seam is solved in the cost matrix by dynamic programming algorithm, and the target image after splicing fusion is output.The present application prevents tearing and overlapping distortion phenomenon in image alignment process.

Description

Technical Field

[0001] This invention relates to the field of image processing, and in particular to a method and system for stitching and fusing graphics and images based on topology analysis. Background Technology

[0002] In practical applications such as industrial image stitching, large-area drone aerial image registration, multi-camera security monitoring fusion, and medical panoramic microscopic image analysis, how to extract reliable features from two images with local overlap and achieve seamless geometric registration has always been a core challenge in the field of computer vision. Existing technologies reveal the following four deep-seated technical bottlenecks: First, false edge contamination caused by complex lighting and smooth backgrounds. Existing technologies typically use a globally fixed grayscale threshold for image binarization; however, images acquired in natural environments often simultaneously contain high-brightness areas exposed to direct sunlight, low-brightness areas with large areas of shadow, and nearly textureless smooth areas such as the sky and white walls. The global threshold easily misclassifies foreground structures in bright areas as background, while in smooth areas, the slight grayscale fluctuations caused by gradual changes in lighting and sensor dark current noise can be incorrectly identified as highly misleading false edges and noise fragments. Once these false edges are mixed in, they will seriously contaminate subsequent structural topology networks. Second, traditional local gradient features are extremely sensitive to large parallax and perspective distortion. Current mainstream local feature extraction algorithms (such as SIFT, SURF, and ORB) all heavily rely on the local pixel gradient distribution of the image in their underlying logic. When the shooting angles of two images change significantly (large parallax), the physical structures in the images will undergo severe affine distortion or perspective distortion, causing traditional feature descriptors to fail and the matching success rate to drop sharply. Third, large parallax mesh deformation optimization suffers from pathological large distortion artifacts due to the lack of hard topological constraints. In large-format geometric mapping caused by large parallax, due to the lack of mathematical constraints on local topological validity, triangular facets are prone to vertex flipping or crossing during numerical optimization iterations, leading to facet orientation reversal or even multiple facets overlapping and penetrating each other. Fourth, seam planning driven by pure color difference is prone to causing visual discontinuity in the structure. In the final image fusion stage, existing seam generation techniques (such as conventional dynamic programming or graph cut algorithms) usually only use the color difference between pixels in the overlapping area as the only pathfinding cost indicator. It is extremely easy to plan a splicing seam that penetrates directly through the middle of the physical subject, producing an extremely abrupt and unnatural visual discontinuity effect. Summary of the Invention

[0003] One of the objectives of this invention is to provide a graphic image stitching and fusion method based on topology analysis to solve the problems of false edge contamination, feature matching failure, mesh folding distortion, and target visual fragmentation in the prior art.

[0004] This invention is achieved through the following technical solution: a graphic image stitching and fusion method based on topological analysis, comprising the following steps: acquiring a first image and a second image; performing binarization partitioning processing on the first image and the second image according to a preset grayscale threshold matrix to generate a first binary image and a second binary image; performing topological feature analysis on the first binary image and the second binary image respectively; generating a deterministic topological feature descriptor by extracting the topological skeleton of the foreground region of the image and calculating the topological characteristic number; generating an initial set of matching point pairs based on the deterministic topological feature descriptor; using the initial set of matching point pairs as data driving terms, dividing the pixel space of the first image into a two-dimensional non-overlapping triangular mesh set, constructing a deformation energy function, injecting Jacobian determinant positive definiteness constraints during the optimization process of the deformation energy function to control the mapping process of the mesh deformation field to satisfy the homeomorphic mapping condition, generating an aligned first image; calculating the overlapping region between the aligned first image and the second image; constructing a cost matrix within the overlapping region based on pixel-level difference values ​​and the repulsion weights of the topological skeleton; solving for the image stitching seam in the cost matrix using a dynamic programming algorithm; and outputting the stitched and fused target image.

[0005] Further, a binarization partitioning process is performed, including: dividing the original image into multiple square blocks, calculating the grayscale mean and variance of each square block; when the variance of a square block is greater than a preset contrast tolerance, setting the binarization threshold of the square block to the grayscale mean; when the variance of a square block is not greater than the preset contrast tolerance, setting the binarization threshold of the square block to a preset background baseline value; comparing the grayscale value of each pixel in the original image with the binarization threshold corresponding to its square block, determining pixels with grayscale values ​​higher than the binarization threshold as foreground structure pixels, and determining pixels with grayscale values ​​not higher than the binarization threshold as background pixels.

[0006] Furthermore, the process of setting the binarization threshold of the square block to switch between the grayscale mean and the background baseline value is implemented through a step switch function. Specifically, it includes: using the difference between the variance of the square block and the contrast tolerance as the input of the step switch function; when the difference is positive, the step switch function outputs one; when the difference is zero or negative, the step switch function outputs zero; multiplying the grayscale mean by the output value of the step switch function, and adding the result of multiplying the background baseline value by one and subtracting the output value of the step switch function, to obtain the dynamic threshold of the square block.

[0007] Further, by extracting the topological skeleton of the foreground region of the image, the process includes: performing a median transformation on the first binary image and the second binary image; generating a distance transformation field by calculating the Euclidean distance from each foreground pixel point inside the image to the nearest background boundary pixel; and connecting the foreground pixels in the distance transformation field that have at least two equidistant nearest neighbor boundary points on the background boundary to form an initial topological skeleton.

[0008] Further, calculating the topological characteristic number includes: abstracting the initial topological skeleton into a one-dimensional simplex complex, wherein the intersection nodes and branch endpoints on the initial topological skeleton are abstracted as vertices, and the skeleton segments connecting adjacent vertices are abstracted as edges; using the Euler-Poincaré formula, calculating the Euler characteristic number represented by the difference between the zeroth-order Betti number and the first-order Betti number based on the difference between the number of vertices and the number of edges of the simplex complex; constructing homology groups based on the boundary operators of the simplex complex, and solving for the zeroth-order Betti number and the first-order Betti number by calculating the rank of each order of homology group, wherein the zeroth-order Betti number represents the number of connected components, and the first-order Betti number represents the number of independent closed holes.

[0009] Furthermore, after forming the initial topological skeleton, the method further includes: traversing all edge branches of the initial topological skeleton, calculating the path integral length of each edge branch, wherein the edge branch is a skeleton line segment segment with the branch endpoint or intersection node as the end, and the path integral length is the cumulative arc length obtained by accumulating the Euclidean distance between all adjacent skeleton pixels along the edge branch; removing edge branches with a path integral length lower than a preset skeleton pruning threshold to obtain a clean topological skeleton; and writing the node coordinates of the clean topological skeleton into the deterministic topological feature descriptor.

[0010] Further, generating a deterministic topological feature descriptor includes: extracting the branch endpoints, intersection nodes, and distance transformation field radius values ​​corresponding to each node from the initial topological skeleton; and concatenating and combining the coordinates of the branch endpoints, the coordinates of the intersection nodes, the distance transformation field radius values ​​of each node, the zeroth Betti number, and the first Betti number to generate the deterministic topological feature descriptor containing structural morphology information and topological invariant information.

[0011] Furthermore, in the process of performing topological feature parsing on the first binary image and the second binary image respectively, the method further includes: when the calculated first-order Betti number and feature endpoint density of the local region are both lower than the preset environmental texture sparsity threshold, a dynamic downscaling analysis mechanism is triggered; the dynamic downscaling analysis mechanism reduces the resolution of the first binary image and the second binary image to one-quarter of the original resolution, expands the receptive field of distance transformation, recalculates the topological skeleton and connected components, and uses the topological consistency of the global macrostructure to replace the locally missing feature alignment points.

[0012] Furthermore, the deformation energy function includes a feature alignment data term and a mesh rigidity preservation regularization term. The feature alignment data term is calculated as follows: for each pair of corresponding points in the initial matching point pair set, the centroid coordinates are interpolated based on the coordinates of the second vertex of the mesh cell where the first image feature point is located in the corresponding point to obtain the mapped position of the first image feature point. The square Euclidean distance between the mapped position and the position of the second image feature point in the corresponding point is calculated, and the square Euclidean distances of all corresponding points are summed. The mesh rigidity preservation regularization term is calculated as follows: for each pair of adjacent vertices in the two-dimensional non-overlapping triangular mesh set, the square Euclidean norm of the deviation between the deformed edge vector and the edge vector before deformation after rotation by the local optimal rotation matrix is ​​calculated. The square Euclidean norms of all adjacent vertex pairs are summed and multiplied by a regularization weight coefficient.

[0013] Furthermore, the positive definiteness constraint of the Jacobian determinant is injected, including: for each mesh cell in the two-dimensional non-overlapping triangular mesh set, calculating the Jacobian matrix corresponding to the local affine mapping of the mesh cell from before deformation to after deformation based on the coordinates of the first vertex before deformation and the coordinates of the second vertex after deformation; calculating the Jacobian determinant of the Jacobian matrix; and forcibly requiring that the Jacobian determinant of each mesh cell in the two-dimensional non-overlapping triangular mesh set is strictly greater than a preset tolerance constant, wherein the tolerance constant is a positive value greater than zero.

[0014] Furthermore, the Jacobian determinant is calculated as follows: the two edge vectors formed by the coordinates of the second vertex of the deformed grid cell are arranged in columns to form a first matrix; the two edge vectors formed by the coordinates of the first vertex of the grid cell before deformation are arranged in columns to form a second matrix; and the inverse of the first matrix and the second matrix are multiplied to obtain the Jacobian matrix of the grid cell.

[0015] Furthermore, a logarithmic obstacle penalty is added to the deformation energy function to enforce the positive definiteness constraint of the Jacobian determinant. The logarithmic obstacle penalty is calculated as follows: for each grid cell in the two-dimensional non-overlapping triangular mesh set, the difference between the Jacobian determinant of the grid cell and the tolerance constant is calculated, and the natural logarithm of the difference is taken; the natural logarithm of all grid cells is summed, multiplied by the obstacle weight coefficient, and the result is negative to obtain the logarithmic obstacle penalty; the deformation energy function and the logarithmic obstacle penalty are summed to form the global optimization objective function.

[0016] Further, generating the aligned first image includes: iteratively updating the coordinates of all second vertices in the two-dimensional non-overlapping triangular mesh set using the Gauss-Newton method; when an iteration causes the Jacobian determinant of any mesh cell to approach the tolerance constant, the logarithmic barrier penalty term tends to positive infinity, and a deterministic gradient backoff mechanism is used to prevent this parameter update, so as to prevent any mesh cell from flipping, overlapping, or tearing; after the iteration converges, the first image is meshed by deformation mapping based on the optimal coordinates of all the second vertices obtained by solving, to generate the aligned first image.

[0017] Further, a cost matrix is ​​constructed based on pixel-level difference values ​​and repulsion weights of the topological skeleton, including: calculating the Euclidean norm of the pixel value difference between the aligned first image and the second image at each pixel position in the overlapping region, generating a difference energy map; introducing the spatial coordinate information of the topological skeleton, calculating the squared value of the Euclidean distance from each pixel point in the overlapping region to the nearest topological skeleton point, summing the squared value with a repulsion weight smoothing constant as the denominator, using a preset positive topological repulsion coefficient as the numerator, calculating the quotient of the numerator and the denominator, and adding the value to the quotient to obtain the repulsion weight; multiplying the Euclidean norm at each pixel position in the difference energy map with the repulsion weight at that position to generate the cost matrix; the preset positive topological repulsion coefficient is used to control the repulsion strength of the topological skeleton to the path of the image stitching seam.

[0018] Further, the image stitching seam is solved in the cost matrix using a dynamic programming algorithm, including: starting from the initial boundary row of the overlapping region, initializing the cumulative energy function value of the row to the cost matrix value of the row; for each pixel position in each subsequent row, adding the cost matrix value of the position to the minimum of the cumulative energy function values ​​of the position and its left and right adjacent positions in the previous row to obtain the cumulative energy function value of the position; when the recursion reaches the last boundary row of the overlapping region, starting from the pixel position with the minimum cumulative energy function value in the row, backtracking from the last boundary row to the first boundary row along the recorded optimal forward direction to obtain the image stitching seam.

[0019] Another aspect of the present invention provides a graphic image stitching and fusion system based on topology analysis, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements any of the graphic image stitching and fusion methods based on topology analysis as described above.

[0020] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0021] 1. This invention breaks through the limitations of global threshold and traditional adaptive threshold. It introduces the local gray-level variance of image blocks as a contrast criterion. Through the hard switching mechanism of nonlinear step switch operator, when the current area is determined to be a smooth background (with minimal variance), the dynamic threshold is forcibly lowered to the background baseline value, eliminating the false edges caused by sensor noise and gradual changes in illumination. This not only separates the extremely pure foreground structure, but also ensures the physical authenticity of the subsequent topological skeleton from the source, cutting off the consumption of computing resources by massive false branches.

[0022] 2. This invention abandons fragile gradient information, introduces algebraic topology theory, extracts a one-dimensional topological skeleton, and uses the Euler-Poincaré formula to calculate the topological characteristic number (multi-order Betti number) that strictly represents the number of connected components and holes. It constructs a deterministic topological feature descriptor that includes geometric properties and topological invariants. As long as the physical objects in the image do not break or stick together, no matter what extreme stretching, twisting, or perspective distortion they undergo, their algebraic topological features maintain strict homeomorphism invariance, improving the matching success rate under large parallax conditions. In addition, to address the problem of mesh folding and flipping caused by large deformation, it innovatively injects the positive definiteness constraint of the Jacobian determinant into the optimization energy function and constructs a topological equivalence penalty term using the logarithmic barrier function. When the Jacobian determinant of any triangular facet approaches the tolerance limit, the penalty term will rise sharply to positive infinity, generating an overwhelming reverse gradient repulsion force, eliminating tearing and overlapping distortion phenomena in the image alignment process. Attached Figure Description

[0023] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0024] Figure 1 This is a flowchart of the method provided in Embodiment 1 of the present invention.

[0025] Figure 2 This is a bar chart comparing the binarization misclassification rates provided in Embodiment 1 of the present invention.

[0026] Figure 3 This is a schematic diagram of the switching behavior curve provided in Embodiment 1 of the present invention.

[0027] Figure 4 This is a comparison diagram of the characteristic stability under different degrees of deformation provided in Embodiment 1 of the present invention.

[0028] Figure 5 This is a schematic diagram comparing the number of pseudo-branches and the feature matching accuracy provided in Embodiment 1 of the present invention.

[0029] Figure 6 This is a schematic diagram of the barrier effect curve provided in Embodiment 1 of the present invention.

[0030] Figure 7 This is a schematic diagram of the distance distribution of the seam path deviating from the skeleton provided in Embodiment 1 of the present invention. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0032] Example 1

[0033] This embodiment discloses a graphic image stitching and fusion method based on topology analysis. Figure 1 The overall method flowchart of this embodiment is shown. As can be seen from the figure, this embodiment includes the following steps:

[0034] Step 1: Obtain the first image and the second image, and perform binarization partitioning on the first image and the second image according to the preset grayscale threshold matrix to generate the first binary image and the second binary image.

[0035] Here, the first image and the second image refer to the two original input images to be stitched and fused. For example, the first image and the second image can be two adjacent aerial photographs with partially overlapping fields of view obtained by drone aerial photography, or two surveillance images with different perspectives covering the same scene taken at the same time by surveillance cameras installed in different locations, or two microscopic images with partial overlap obtained by scanning the same tissue section with a microscope field by field in the field of medical imaging.

[0036] A grayscale threshold matrix is ​​a pixel-by-pixel defined threshold field with the same spatial resolution as the image. Each element in this threshold field corresponds to a binary decision threshold at the same spatial location in the original image. This grayscale threshold matrix is ​​not a globally uniform fixed constant, but rather a spatially adaptive threshold distribution dynamically calculated based on the local illumination characteristics and texture contrast features of different spatial regions within the image.

[0037] Binarization partitioning refers to using a grayscale threshold matrix to perform pixel-by-pixel threshold comparison and decision operations on the original grayscale image. The grayscale value of each pixel in the original image is compared with the local dynamic threshold at its corresponding spatial location. If the grayscale value is higher than the threshold, the pixel is determined to be a foreground structure pixel and assigned a value of 1; otherwise, the pixel is determined to be a background pixel and assigned a value of 0. In this way, the original image containing complex grayscale levels is transformed into a binary image containing only foreground and background values.

[0038] The first binary image and the second binary image refer to the clean foreground-background separated images obtained after performing the above-mentioned binarization partitioning process on the first image and the second image, respectively, in which only the truly effective foreground structure information in the original image is retained.

[0039] In this embodiment, the methods for acquiring the first image and the second image include, but are not limited to: acquiring images sequentially by means of an aerial camera mounted on a drone during continuous flight according to a preset overlap rate; acquiring images synchronously by means of a multi-camera fixed monitoring system; acquiring images by means of a stepper motor driven by a microscope stage to scan each field of view; or acquiring images directly from a pre-stored image database.

[0040] In this embodiment, in actual industrial-grade image stitching and registration applications, the input images are often acquired under extremely complex and uncontrollable natural environments. For example, when a drone performs aerial photography missions over a city, the images captured by its lens simultaneously include building facades directly exposed to sunlight (high-brightness areas), large areas of shadow cast by the back of buildings (low-brightness areas), and large smooth areas such as the sky and white walls that contain almost no texture details. In these large smooth areas, due to unavoidable physical factors such as gradual changes in illumination and dark current noise from sensors, pixel values ​​are not strictly constant but exhibit weak and slow grayscale fluctuations. If a globally fixed grayscale threshold is directly used to binarize the entire image, a large number of foreground pixels may be incorrectly identified as background in high-brightness areas, while background noise may be incorrectly elevated to foreground structures in low-brightness areas. More seriously, in smooth areas such as the sky or white walls, the weak grayscale fluctuations caused by gradual changes in illumination will be incorrectly identified as valid edges or textures by the global threshold, resulting in a large number of false edges and noise fragments in the binarization result. Once this false structural information is introduced into subsequent topological feature extraction and skeleton analysis, it will directly lead to a large number of meaningless noise branches in the topological network, severely polluting the purity of the topological features and ultimately causing a sharp decline in the accuracy and robustness of the entire registration and stitching process. Therefore, before performing any subsequent feature extraction and topological analysis, the first and crucial step is to design an intelligent binarization mechanism that can adapt to the differences in illumination characteristics in different regions within the image, eliminating the interference of noise on subsequent processing links from the source.

[0041] In this embodiment, to overcome the inherent limitations of the global threshold, the grayscale threshold matrix preset in step S1 is a spatially adaptive threshold dynamically generated based on the local variance of image blocks, specifically including:

[0042] The original image is divided into multiple square blocks, and the grayscale mean and variance of each block are calculated. When the variance of a block is greater than the preset contrast tolerance, the binarization threshold of the block is set to the mean. When the variance of a block is less than the preset contrast tolerance, the threshold of the block is forcibly shifted to the background baseline value, thereby filtering out the illumination gradient interference in large smooth areas.

[0043] Among them, a square block refers to dividing the original image into several non-overlapping local square regions in space according to a preset square side length.

[0044] The grayscale mean is a scalar value obtained by taking the arithmetic mean of the grayscale values ​​of all pixels contained in each square block. This mean reflects the overall brightness level of the block under the current lighting conditions.

[0045] Variance is a scalar value obtained by taking the arithmetic mean of the squared deviations of the gray values ​​of all pixels within a square block from the mean gray value of that block. Variance is a key statistic for measuring the dispersion of pixel values ​​and the strength of contrast within a region: if the gray variance of a region is large, it indicates that there are significant changes in brightness and darkness within the region, and it is likely to contain effective foreground structure edge information; conversely, if the gray variance of a region is extremely small, it indicates that the region is a nearly uniform and smooth area (such as the sky or a white wall), and the slight gray fluctuations within it are most likely due to gradual changes in lighting or sensor noise, and should not be considered as effective foreground structure.

[0046] The preset contrast tolerance refers to a pre-defined variance critical criterion constant used to distinguish between effective texture regions and smooth background regions. When the variance of a region is lower than this constant, the model determines that the region is a smooth background region that does not contain effective foreground information. The value of this constant can be calibrated according to the noise level of the imaging device and the intensity of gradual illumination changes in the specific application scenario.

[0047] The background baseline value is a pre-set grayscale constant, usually set to a sufficiently low value or an estimated value based on the global background grayscale level. When a region is identified as a smooth background region, the dynamic threshold of that region will be forced to be set to this background baseline value, thereby ensuring that all pixels in that region are classified as background when compared with this threshold, thus completely eliminating false edges in smooth regions.

[0048] Specifically, the original input image It is spatially divided into a set of non-overlapping squares. For any square region arbitrary pixel coordinates within The binarization discrimination result is determined by comparing the gray value of the pixel with a dynamic threshold specific to the region.

[0049] For example, in this embodiment, the binarization discriminant function can be expressed as follows:

[0050]

[0051] in, pixel coordinates The binarized output at the point is such that when its value is 1, it indicates that the pixel is determined to be a foreground structure pixel, and when its value is 0, it indicates that the pixel is determined to be a background pixel. The original input image in coordinates The grayscale value at that location; For pixels Belonging to the Each square area The corresponding dynamic threshold is determined adaptively by the local statistical characteristics of the pixels within the square region, rather than being globally uniform and fixed.

[0052] The dynamic threshold matrix in the above discriminant function It is the core of the entire adaptive binarization model. The calculation of this threshold matrix is ​​determined by the local statistics (mean and variance) of the region and a nonlinear step operator. Its design goal is to achieve automatic masking of false edges in smooth regions.

[0053] For example, in this embodiment, the dynamic threshold matrix can be calculated using the following formula:

[0054]

[0055] in, For the first Each square area The local mean of all pixel grayscale values ​​within the region reflects the overall brightness level of the area under the current lighting conditions. In areas containing effective foreground structures, it will serve as a reference grayscale line for judging the foreground and background. For the first Each square area The local variance of all pixel grayscale values ​​within a region is a key indicator for measuring the degree of grayscale dispersion and contrast strength of pixels within that region. The larger the variance value, the more significant the alternation of light and dark in that region, and the more likely it is to contain effective foreground structure boundary information. Conversely, a very small variance value means that the region is a nearly uniform and smooth area. The Heaviside step function is a classic nonlinear switching operator. It is defined as follows: when the value inside the parentheses is positive, the function output is always equal to 1; when the value inside the parentheses is zero or negative, the function output is always equal to 0. In this model, the step function plays the core role of logical branch judgment, and decides to switch the threshold calculation mode based on the relationship between the local variance and the tolerance threshold. This is the preset contrast tolerance constant; This is the background baseline value. Figure 2 A bar chart comparing the binarized misclassification rates of the spatial adaptive threshold and the global fixed threshold in this embodiment is shown. Figure 2 The X-axis represents the image block index, and the Y-axis represents the foreground false positive rate (%) for that block. The global threshold significantly misses foreground detections in high-brightness areas and significantly misdetects background noise in low-brightness areas, while the adaptive threshold maintains an extremely low false positive rate across all regions. Figure 2 As can be seen, the block dynamic threshold model has the ability to adaptively compensate for illumination non-uniformity, eliminating false edge interference from the source. Figure 3This embodiment illustrates the switching behavior curves of the dynamic threshold under different block variances. Figure 3 A step function-driven intelligent switching mechanism was demonstrated, verifying the complete shielding effect on pseudo-edges in smooth regions.

[0056] It is understandable that in the above formula, This constitutes the intelligent switching mechanism of the entire dynamic threshold model. When a certain cube area is located in a high-contrast area with significant texture structure, such as a building facade or road edge, the gray-level variance of that area... Much larger than the tolerance constant Step function The output is 1, at which point the dynamic threshold is... In other words, the threshold is dominated by the local mean of the region, and foreground-background separation is performed according to normal adaptive logic. However, when a certain square region is located in a large, smooth area such as the sky or a white wall, the grayscale variance of that region... Much smaller than the tolerance constant The step function output is 0, at which point the dynamic threshold... In other words, the threshold is forcibly switched to the background baseline value, suppressing all subtle gray-level fluctuations in the region caused by gradual illumination changes as background. Through this hard-switching mechanism based on a step function, the model fundamentally achieves automatic discrimination and differential processing between meaningful textured regions and meaningless smooth regions. This ensures that the final binarized image retains only truly effective foreground structural information, providing a crucial prerequisite for extracting a clean topological skeleton and effectively avoiding a large number of false branches caused by noise in the subsequent topological network.

[0057] Step 2: Perform topological feature parsing on the first binary image and the second binary image respectively. Generate deterministic topological feature descriptors by extracting the topological skeleton of the image region and calculating the topological characteristic number. Then, compare the topological feature descriptors to generate an initial set of matching point pairs.

[0058] Topological feature analysis refers to the process of extracting strictly invariant topological invariants from the foreground structure of a binary image using theoretical tools of algebraic topology. This analytical process focuses not on the specific geometric shape or size of the object, but rather on its topological properties that remain unchanged under continuous transformations—such as the number of independent connected components or closed holes within a structure.

[0059] A topological skeleton is a one-dimensional linear central axis structure extracted from a two-dimensional foreground region through a central axis transformation. This skeleton faithfully preserves the topological connectivity and overall morphological information of the original two-dimensional foreground region in a one-dimensional linear representation.

[0060] Topological characteristic numbers are algebraic invariants that describe the essential characteristics of topological structures, calculated using the Euler-Poincaré formula. Specifically, they include the 0th-order Betti number. and the first-order Betty number .

[0061] A deterministic topological feature descriptor is a multidimensional feature vector that combines the geometric properties of the topological skeleton (branch endpoints, intersection nodes and their corresponding distance transformation field radii) with topological invariants (0th-order Betti number and 1st-order Betti number). This feature vector simultaneously encodes the structural morphology information of the image foreground and the topological invariant information.

[0062] The initial matching point pair set refers to a set of topologically consistent feature correspondences established between two images by comparing the deterministic topological feature descriptors generated by the first binary image and the second binary image.

[0063] In this embodiment, the most widely used local feature extraction methods (such as SIFT, SURF, ORB, etc.) in traditional image registration and matching techniques essentially rely on pixel-level gradient information in the image to detect and describe key points. However, these gradient-based feature descriptors have a fundamental theoretical flaw: they are extremely sensitive to geometric deformations such as affine transformations and large parallax perspective transformations. When the shooting angle changes significantly, structural features that appear as perfect circles in one image may be distorted into extremely flat ellipses in another image due to perspective effects. In these cases, the values ​​of local gradient-based feature descriptors will change drastically, leading to a sharp drop in the success rate of feature matching or even complete failure. To fundamentally overcome this bottleneck, this method introduces core theoretical tools from algebraic topology to extract structural features of the image from a higher-dimensional mathematical abstraction level. The Betti number combination extracted by algebraic topology... As a topological invariant, it is strictly invariant under homeomorphic transformations and naturally possesses perfect immunity to geometric deformations.

[0064] In this embodiment, step S2, which involves extracting the topological skeleton of the image region and calculating the topological characteristic number to generate a deterministic topological feature descriptor, further includes:

[0065] S2.1: Perform median transformation on the first binary image and the second binary image, generate a distance transformation field by calculating the Euclidean distance from each non-zero pixel point inside the image to the nearest background boundary, and connect the local maxima points in the distance transformation field to form an initial topological skeleton.

[0066] Among them, the median transformation is a classic morphological analysis method. Its core modeling idea is to calculate the Euclidean distance from each foreground pixel in the binary image to the nearest background boundary pixel, thereby constructing a continuous distance field over the entire foreground region.

[0067] The distance transform field is a scalar field formed by calculating the Euclidean distance from each foreground pixel to the nearest background boundary over the entire foreground region of a binary image. In this distance field, the pixel located at the center of the foreground region has the largest distance value because it is farther from the boundaries in all directions; while the closer the pixel is to the foreground boundary, the smaller its distance value.

[0068] The initial topological skeleton refers to the ridges or singularity set in the distance transform field, that is, the set of foreground pixels that have two or more equidistant nearest neighbors on the background boundary. The points on these ridges have a significant geometric property: they have at least two equidistant nearest neighbors on the background boundary.

[0069] Specifically, the binary domain obtained after step 1 First, determine the non-zero foreground subset consisting of all pixels with a value of 1. and the set of boundary pixels of the foreground subset. Then, for the foreground subset Each pixel within Calculate the distance from the pixel to the boundary set. Calculate the Euclidean distances of all boundary pixels in the range, and take the minimum value among them as the distance transformation value of that point.

[0070] For example, in this embodiment, the Euclidean distance field can be calculated using the following formula:

[0071]

[0072] in, Foreground pixels The distance transformation value at a given location represents the Euclidean distance between the pixel and the nearest background boundary pixel. A larger value indicates that the pixel is located deeper inside the foreground region, while a smaller value indicates that the pixel is closer to the edge contour of the foreground. The physical meaning of the bounding operator here is the boundary set. Search and return all boundary pixels with the same value as the current foreground pixel. The value with the smallest distance; For traversing the boundary set The coordinate variables used for each boundary pixel; This is a two-dimensional Euclidean norm operator used to calculate the straight-line distance between two pixel coordinates.

[0073] After obtaining the above distance transformation field Then, the initial topological skeleton of the foreground region can be extracted based on this.

[0074] For example, in this embodiment, the mathematical definition of the topological skeleton can be shown as follows:

[0075]

[0076] in, The extracted initial topological skeleton set is a one-dimensional central axis representation of the original two-dimensional foreground region, faithfully preserving the topological connectivity and overall morphological structure of the foreground region; Foreground subset Any pixel in the array that is being examined; and Boundary set Two distinct boundary pixels; condition Represents pixels There exist at least two distinct nearest neighbor boundary points on the background boundary. and And the distance to both boundary points is exactly equal to the distance transformation value of that point. Points that satisfy this condition are located on the ridge line of the distance transformation field and are part of the central axis of the foreground region.

[0077] Understandably, the essence of the above-described topological skeleton extraction process is to compress a two-dimensional planar foreground region into a one-dimensional linear skeleton structure, while completely preserving the topological connections of the original foreground region. For example, a complete annular foreground region, after a mid-axis transformation, will yield a closed annular skeleton line, which faithfully reflects the topological feature that the original annulus has a hole in the middle; while two separate foreground regions will yield two unconnected skeleton line segments, reflecting the topological feature that there are two independent connected branches.

[0078] S2.2 Constructing a simple complex data structure involves triangulation of the initial topological skeleton, and calculating the 0th order Betti number within the region using the Euler-Poincaré formula. and the first-order Betty number ,in Indicates the number of connected components. This indicates the number of individual holes.

[0079] Simplex is a fundamental data structure in algebraic topology, which performs algebraic analysis by discretizing geometric objects into a combination of topological structures consisting of vertices and edges.

[0080] Triangulation refers to the process of abstracting an initial topological skeleton into a one-dimensional simplex complex, where the intersections and endpoints of the skeleton lines are abstracted as vertices, and the skeleton line segments connecting adjacent vertices are abstracted as edges.

[0081] The Euler-Poincaré formula is a classic formula in algebraic topology that establishes a bridge between geometric quantities (number of vertices, number of edges) and topological invariants (Betti numbers).

[0082] 0th order Betty number Strictly defined algebraically as the rank of a zero-dimensional homology group, its physical meaning is the number of independent connected components in a binary foreground of an image. For example, when there are three separate buildings in the image... .

[0083] First-order Betty number Strictly defined algebraically as the rank of a one-dimensional homology group, its physical meaning is the number of closed holes in the binary foreground of an image. For example, a hole enclosed by a window frame or a hole in the middle of a ring will contribute a rank. count.

[0084] Specifically, the skeleton Abstracted as a one-dimensional simple complex This simple complex is composed of a set of vertices. (Intersections and endpoints on the skeleton line) and edge sets (The skeleton segments connecting adjacent vertices) constitute the complex. By calculating the Euler-Poincaré characteristic of this simplex, an algebraic relationship is established between the number of vertices, the number of edges, and the Betti number.

[0085] For example, in this embodiment, the Euler-Poincaré characteristic formula can be expressed as follows:

[0086]

[0087] in, For simple complex The Euler-Poincaré characteristic is an integer-valued topological invariant that represents the most fundamental algebraic characteristic of the topological structure. For the vertex set of a simple complex The cardinality, that is, the total number of intersections and endpoints in the skeleton structure; For edge sets in simple complex The cardinality is the total number of skeleton segments connecting adjacent vertices; It is the 0th order Betty number. It is a first-order Betty number.

[0088] To solve separately and The specific value needs to be further introduced by introducing the concept of boundary operators from the theory of simple homology, and the Betti number of each order is strictly defined by calculating the rank of the homology group.

[0089] For example, in this embodiment, the Betti numbers of each order can be calculated by the following formula:

[0090]

[0091] in, For the first in simple homology theory The order boundary operator is a... 3D simplex mapping to its The linear mapping of the dimensional boundary can be represented as an integer matrix during the computation process. This matrix encodes the correlation between the simplexes of different dimensions in the simplex complex. For the first The kernel space of the order boundary operator, i.e., all boundaries are zero. A collection of dimensional chains; This is the image space of the second-order boundary operator; Operators for determining the dimension of a vector space.

[0092] Figure 4 This diagram shows a comparison of the Betti number and the traditional gradient feature (SIFT) under different degrees of deformation in this embodiment. Figure 4 The X-axis represents the perspective deformation intensity (tilt angle from 0° to 80°), and the Y-axis represents the normalized consistency score (0~1) of the feature descriptor. The consistency of SIFT feature descriptors decreases sharply with increasing deformation angle, dropping to almost zero after exceeding 40°; while the Betty number combination... The fact that it remains essentially constant across the entire angular range demonstrates the strict invariance of algebraic topological invariants under homeomorphic transformations, proving that our scheme has a fundamental theoretical advantage over traditional gradient features in scenarios with large parallax.

[0093] Understandably, Betty counts and As algebraic topological invariants, they possess a crucial theoretical property in this application: they are strictly invariant under homeomorphic transformations. This means that no matter how severely objects in an image are tilted, stretched, compressed, or perspective-distorted (e.g., a perfect circle is deformed into an extremely flattened ellipse by extreme perspective effects), as long as the object does not physically break or merge, its corresponding Betti number combination remains unchanged. It will remain absolutely unchanged.

[0094] S2.3 Extract the branch endpoints, intersection nodes and their corresponding distance transformation field radius values ​​of the initial topological skeleton, and concatenate them with the 0th order Betti number and the 1st order Betti number to generate a deterministic topological feature descriptor that includes structural morphology and topological invariant dimensions.

[0095] In this context, a branch endpoint refers to a terminal node in the initial topological skeleton that is connected to only one skeleton line segment, i.e., a vertex with a degree of 1.

[0096] A cross node is a branching node in the initial topology skeleton that connects three or more skeleton segments, i.e., a vertex with a degree greater than or equal to 3.

[0097] The radius value of the distance transformation field refers to the distance transformation field between the branch endpoints and the intersection nodes. The corresponding distance value reflects the local width of the foreground region represented by the skeleton node at that location.

[0098] Understandably, deterministic topological feature descriptors encode both the local geometric properties of the skeleton (endpoint coordinates, node coordinates, radius values, and other structural morphological information) and global topological invariants (…). , This achieves the ability to represent both local details and global topological constraints in a unified feature vector.

[0099] In this embodiment, trimming determination logic may also be included after S2.3:

[0100] S2.4. Traverse all edge branches of the initial topological skeleton and calculate the path integral length of each edge branch.

[0101] Among them, edge branches refer to skeleton line segments in the initial topology skeleton with the branch endpoint as one end and the intersection node as the other end (or with both endpoints as the two ends).

[0102] The path integral length refers to the cumulative arc length along the skeleton line of a certain edge branch from the starting node to the ending node. This value is obtained by summing the Euclidean distances between all adjacent skeleton pixels on the branch segment by segment.

[0103] S2.5 Remove edge branches whose path integral length is lower than the preset skeleton trimming threshold to eliminate pseudo skeleton branches caused by high-frequency noise at the image edges, and obtain a pure topological skeleton. Then, write the node coordinates of the pure topological skeleton into the topological feature descriptor.

[0104] The preset skeleton trimming threshold is a pre-defined length threshold used to distinguish between effective skeleton branches generated by the real foreground structure and false short branches generated by high-frequency noise at the image edges.

[0105] Pseudo-skeleton branches refer to extremely short, physically meaningless skeleton spikes generated during the median transformation due to jagged high-frequency noise at the edges of the binarized image.

[0106] A pure topological skeleton refers to a simplified skeleton structure that, after the above pruning process, retains only the effective skeleton branches whose path integral length is greater than the pruning threshold.

[0107] Understandably, the purpose of the trimming decision logic is to improve the signal-to-noise ratio of the topological skeleton. In actual binarized images, the edge contours of the foreground region are often not ideally smooth curves, but rather contain tiny jagged edges caused by pixelation discretization effects and high-frequency noise. These tiny jagged edges generate a large number of extremely short spurious skeleton branches (i.e., burrs) during the median transformation. If these burrs are not cleaned up, they will be mixed into the topological feature descriptor as false endpoints and branches, severely interfering with the accuracy of subsequent feature matching. By setting a lower threshold for the path integral length and removing all short edge branches below this threshold, these noise burrs can be effectively removed, resulting in a pure topological skeleton that only reflects the true macroscopic structure of the foreground. Figure 5 This diagram illustrates a comparison of the number of pseudo-branches and feature matching accuracy before and after skeleton pruning in this embodiment. Figure 5 The chart is a dual Y-axis bar chart. The X-axis represents different noise levels (low / medium / high / extremely high). The left Y-axis represents the number of pseudo branches in the skeleton (bars), and the right Y-axis represents the feature matching accuracy (%, line). The number of pseudo branches in the unpruned skeleton increases sharply with the increase in noise, causing the matching accuracy to plummet. After pruning, the number of pseudo branches is almost zero, and the matching accuracy remains at a high level.

[0108] In this embodiment, step 2 may further include adaptation extensions for complex scenarios:

[0109] When the calculated local region's first-order Betti number When both the density of feature endpoints and the density of environmental textures are lower than the preset environmental texture sparsity threshold, a dynamic downscaling analysis mechanism is triggered.

[0110] Among them, the environmental texture sparsity threshold refers to a pre-set combined criterion used to determine whether the topological features of the current image region are too sparse to provide enough local feature alignment points.

[0111] The dynamic downscaling analysis mechanism refers to the system reducing the analytical resolution of the first binary image and the second binary image to one-quarter of the original resolution when it is determined that the topological features at the current resolution are too sparse. This expands the receptive field of the local Euclidean distance transformation, recalculates the topological skeleton and connected components of the global contour, and uses the topological consistency of the global macrostructure to replace the locally missing feature alignment points.

[0112] Understandably, reducing the resolution to one-quarter is equivalent to performing a 2x spatial downsampling on the binary image. This operation increases the effective physical area covered by each pixel, thus doubling the receptive field of the distance transform. In the downsampled low-resolution space, large-scale macroscopic contour structures (such as the overall outline of buildings and the boundaries of large features) that were previously impossible to capture effectively at high resolution due to their small spatial extent will be fully presented in the topological skeleton. By recalculating the topological skeleton and Betti number at low resolution, the system can establish a global correspondence between the two images using the topological consistency of these large-scale macroscopic structures, thereby maintaining the effectiveness of registration constraints even in scenes with severely lacking local detail features.

[0113] Step 3: Using the initial set of matching point pairs as the data driver, construct a mesh deformation field from the first image coordinate system to the second image coordinate system. In the optimization process of the mesh deformation field, inject Jacobian determinant positive definiteness constraints to control the mapping process of the mesh deformation field to satisfy the homeomorphic mapping condition and generate the aligned first image.

[0114] Among them, the mesh deformation field refers to the continuous spatial mapping relationship of the source image from its original geometric state to the target alignment state by adjusting the coordinate positions of each vertex in the triangular mesh after parameterizing the two-dimensional plane domain where the first image is located into a triangular mesh structure.

[0115] The Jacobian determinant positive definiteness constraint is a mathematical hard constraint that requires the determinant value of the Jacobian matrix in every local region of the mesh deformation field to always be strictly positive.

[0116] The homeomorphic mapping condition refers to the topological equivalence requirement that a deformation mapping must satisfy: the mapping must be continuous and invertible, and its inverse mapping must also be continuous, thereby ensuring that no local flipping, overlapping, tearing, or dimensional collapse occurs during the deformation process.

[0117] In this embodiment, in the engineering practice of image registration, when facing image registration tasks under conditions of large parallax (e.g., photographing the same building from two far apart camera positions), the geometric differences between the source and target images are often extremely large. This requires that certain local areas in the mesh must undergo very drastic deformation to achieve alignment. Under such large deformation conditions, when traditional registration energy functions are used for numerical optimization, the relative positions of the vertices of some triangular patches may flip or cross during the iterative update process, leading to a reversal of the geometric orientation of the triangular patches, and even severe degradation such as multiple triangular patches overlapping and penetrating each other. This mesh folding and flipping phenomenon visually manifests as severe distortion artifacts in the registered image, such as buildings appearing to be broken in half and texture areas overlapping and wrinkling. The root cause of this problem is that the traditional registration energy function is a purely quadratic function, which does not impose any mathematical constraints on the topological validity of the mesh during the deformation process.

[0118] To fundamentally solve the above problems, this method introduces an additional topological constraint based on the positive definiteness of the Jacobian determinant in differential geometry, on the basis of the traditional registration energy function. By using the interior point method logarithmic barrier function technique in numerical optimization theory, this hard constraint is transformed into a continuously differentiable penalty term with an absolutely rigid mathematical barrier effect. This ensures that every triangular facet of the mesh always maintains a positive orientation throughout the entire optimization iteration process and will never fold or flip, thus strictly guaranteeing the homeomorphism of the deformation mapping mathematically.

[0119] In this embodiment, step 3, which uses the initial set of matched point pairs as data-driven terms to construct the mesh deformation field from the first image coordinate system to the second image coordinate system, may further include:

[0120] S3.1. The pixel space of the first image is uniformly divided into a set of two-dimensional non-overlapping triangular meshes composed of multiple first vertices.

[0121] Here, the first vertex refers to the original two-dimensional coordinate position of the corner vertex of each triangular facet in the triangular mesh set before deformation.

[0122] A two-dimensional non-overlapping triangular mesh set refers to a planar partitioning structure that covers the entire pixel space of a first image and is composed of multiple non-overlapping triangular facets. Each triangular facet in this mesh set shares a boundary with its adjacent facets, and there are no overlaps or gaps between facets.

[0123] S3.2 Based on the coordinate correspondence of the initial matching point set, construct a deformation energy function to calculate the coordinates of the second vertex corresponding to the deformation of the two-dimensional non-overlapping triangular mesh set. The deformation energy function includes a feature alignment data term and a mesh rigidity preservation regularization term.

[0124] Here, the second vertex refers to the two-dimensional coordinates of the target object after deformation of all vertices in the mesh, which is the decision variable of the entire optimization problem.

[0125] The deformation energy function is a scalar objective function. The smaller the value, the closer the current mesh deformation state is to the desired alignment target.

[0126] The feature alignment data term refers to the part of the deformation energy function that measures the geometric deviation between matched feature point pairs. This term drives the mesh to deform in the direction that reduces the matching error.

[0127] The mesh rigidity preservation regularization term refers to the part of the deformation energy function that constrains the relative positional relationship between adjacent vertices to be as close as possible to a local rigid rotation transformation, in order to prevent the mesh from undergoing excessively free unconstrained deformation.

[0128] For example, in this embodiment, the triangular mesh of the first image domain is... The set of vertex vectors in is The deformed target coordinates are as follows The deformation energy function can be expressed as follows:

[0129]

[0130] in, This is the deformation energy function; the smaller the value, the closer the current mesh deformation state is to the desired alignment target. For all in the grid A matrix consisting of the original two-dimensional coordinates of each first vertex, each row... Representing the Initial coordinates of each vertex ; Let be the matrix consisting of the two-dimensional coordinates of the target object at all the second vertices in the grid. This matrix represents the decision variables for the entire optimization problem. This is the transpose symbol for a matrix; Represents the initial set of matching point pairs A pair of corresponding points in, where The location of the feature point in the first image. For the second image and The corresponding feature point locations; To be determined by the current mesh deformation state The defined mapping function is based on feature points. The vertex coordinates of the triangular facet are calculated using barycentric coordinate interpolation. At the new position after deformation This is the squared Euclidean norm of the alignment error of the feature point pair, which is the specific implementation of the feature alignment data item; The regularization weight coefficient is used to balance the relative importance between the feature alignment data term and the mesh rigidity preservation regularization term. The larger the value, the more rigid the deformation tends to be, and the smaller the value, the greater the degree of flexible deformation is allowed. This indicates traversing all adjacent vertex pairs in the grid; For defined on the edge The local optimal rotation matrix on the edge describes the ideal rigid rotation angle that the edge should undergo within the local region. This measure is the degree of deviation between the edge vectors of adjacent vertices after deformation and the edge vectors after ideal rigid rotation. This term is the specific implementation of the mesh rigidity preservation regularization term.

[0131] It is understandable that the above deformation energy function Although it can effectively drive the mesh to optimize towards data alignment and local rigidity, it is itself an unconstrained quadratic energy function. When performing numerical optimization under large parallax conditions, there is no mathematical mechanism to prevent the relative position flip of vertex coordinates of certain triangular patches during the iteration process.

[0132] S3.3 Establish a set of vertex coordinate transformation equations to parameterize the mesh deformation field into unknown coordinate vectors of all second vertices.

[0133] The vertex coordinate transformation equations refer to the set of derivatives of the deformation energy function with respect to the coordinates of all second vertices. This set of equations transforms the optimization problem of the mesh deformation field into a numerical optimization problem with the coordinates of all second vertices as unknown variables.

[0134] Understandably, this is achieved by parameterizing the mesh deformation field as a vector containing the unknown coordinates of all second vertices. The entire registration problem is uniformly formulated as a nonlinear optimization problem concerning the vector, which facilitates the subsequent joint solution by introducing the positive definiteness constraint of the Jacobian determinant and the logarithmic obstacle penalty.

[0135] In this embodiment, step S3, which involves injecting Jacobian determinant positive definiteness constraints during the optimization process of the mesh deformation field to control the mapping process of the mesh deformation field to satisfy the homeomorphic mapping condition, specifically includes:

[0136] S3.4 For any mesh element in a two-dimensional non-overlapping triangular mesh set, let the local coordinates before deformation be... The transformed continuous mapping function is Calculate the Jacobian matrix corresponding to the deformation. .

[0137] In this context, a mesh cell refers to any triangular facet in a set of two-dimensional non-overlapping triangular meshes.

[0138] Local coordinates This refers to the coordinates of the triangular facet in the reference coordinate space before deformation.

[0139] Continuous mapping function This refers to the local affine mapping relationship of the triangular facet from the reference space before deformation to the target space after deformation.

[0140] Jacobian matrix It is The matrix completely describes the local linear transformation characteristics that the triangular facet undergoes during deformation, including rotation, scaling, and shearing.

[0141] Specifically, for grids Any triangular facet in Let the coordinates of its three vertices before deformation be... The coordinates of the three vertices after deformation are The geometric transformation of the triangular facet from before to after deformation can be described by a local affine mapping, and the Jacobian matrix of this local affine mapping can be explicitly calculated by the algebraic relationship between the vertex coordinates before and after deformation.

[0142] Exemplarily, in this embodiment, a triangular facet The Jacobian matrix of a local affine mapping can be calculated using the following formula:

[0143]

[0144] in, Triangular facet Local affine mapping Jacobian matrix; It is formed by arranging the two side vectors of the deformed triangle in columns. matrix; It is formed by arranging the two side vectors of the original triangle in columns. The inverse of the matrix is ​​a fixed, known quantity that only needs to be calculated once before optimization begins.

[0145] S3.5 Solving for the determinant of the Jacobian matrix .

[0146] The Jacobian determinant is the signed ratio of the area of ​​the deformed triangle to the area of ​​the triangle before deformation. A positive value indicates that the triangle maintained its original vertex orientation during deformation (e.g., always counterclockwise), meaning no flipping occurred. A negative value indicates that the vertex orientation of the triangle was reversed, meaning the triangle folded or flipped during deformation. A zero value indicates that the triangle degenerates to zero area (three vertices collinear), meaning dimensional collapse occurred.

[0147] It is understandable that the Jacobian determinant There is a strict mathematical equivalence between the positive definiteness of a mapping and the homeomorphism of a mesh deformation mapping. According to the inverse function theorem in differential geometry, for a continuously differentiable mapping, if its Jacobian matrix is ​​non-singular at every point in the domain (i.e., the determinant is not zero), and the sign of the determinant remains consistent throughout the domain (i.e., always positive or always negative), then the mapping is locally a homeomorphism. Therefore, requiring that the Jacobian determinant of all triangular faces be strictly positive is mathematically equivalent to guaranteeing the homeomorphism and topological validity of the entire mesh deformation mapping.

[0148] S3.6. Transform the global deformation optimization problem into an extremum problem with inequality constraints, forcing the mathematical constraints to be satisfied for any mesh element: ,in It is a preset tolerance constant greater than zero, thereby ensuring that non-topological equivalence phenomena such as local flipping, overlapping or tearing are strictly prohibited in adjacent mesh areas after deformation.

[0149] Among them, the preset tolerance constant It is a pre-defined positive baseline constant that ensures that even under the most extreme deformation conditions, the Jacobian determinant of each triangular facet will not degenerate to zero, but will always remain at a constant value. Within the above safe positive range, the occurrence of zero-area degradation is completely eliminated from a numerical perspective.

[0150] Non-topological equivalence phenomena refer to pathological distortions that disrupt the original topology during mesh deformation, such as local flipping (triangular facets facing in reverse), overlapping (multiple facets penetrating each other), or tearing (mesh connectivity breaking).

[0151] In this embodiment, to strictly enforce the constraints, a topological equivalence penalty term is introduced using a logarithmic barrier function during the solution process:

[0152] S3.7, The form of the logarithmic barrier penalty added to the deformation energy function is as follows: ,in For the first Jacobian matrix of grid cells This is the penalty weighting coefficient.

[0153] Among them, the logarithmic barrier penalty is a continuously differentiable penalty term constructed based on the classic interior-point method in numerical optimization—the logarithmic barrier function. The basic principle of the logarithmic barrier function is: the logarithmic function... exist Approaching zero, it tends towards negative infinity. Utilizing this mathematical property, a penalty term can be constructed such that when the Jacobian determinant of a triangular facet approaches the tolerance threshold (i.e., approaches folding), the value of this penalty term rapidly approaches positive infinity, thus forming an infinitely high energy barrier in the gradient field of the objective function, preventing the optimization algorithm from pushing any triangular facet into a folded state.

[0154] Penalty weighting coefficient Also known as the barrier parameter, it controls the relative importance of the logarithmic barrier term in the total energy function. In practical solutions, a gradually decreasing annealing strategy is typically used to adjust this parameter—setting a larger value in the early stages of optimization. The initial value is set to ensure that all faces are far from the fold boundary, thus obtaining a globally safe initial solution. This value is then gradually reduced in the later stages of optimization. This value allows for more precise positional adjustments to the face within a safe range, thereby improving alignment accuracy.

[0155] For example, in this embodiment, the topological equivalence penalty can be calculated using the following formula:

[0156]

[0157] in, This is the total energy value of the topological equivalence penalty. The purpose of this term is to completely prevent any triangular facets from folding or flipping during mesh deformation from a mathematical perspective. This refers to the barrier weight coefficient; Triangular facet The determinant value of the Jacobian matrix is ​​physically represented by the signed area scaling factor of the patch. This is a preset tolerance constant; For the natural logarithm function, when the value inside the parentheses... When approaching zero, It tends toward negative infinity, the preceding negative sign Make It tends towards positive infinity, thus creating an infinitely high barrier in the energy field; summation symbol This means that the above logarithmic barrier constraint is applied independently to each triangular facet in the mesh, ensuring that the topological validity of each local region is protected without omission.

[0158] The above deformation energy function and topological equivalence penalty are integrated into a global nonlinear homeomorphic optimization model.

[0159] For example, in this embodiment, the overall objective function of the global nonlinear homeomorphic optimization model can be expressed as follows:

[0160]

[0161]

[0162] in, The optimal solution to the global optimization problem is the optimal post-deformation coordinates of all mesh vertices under the condition of satisfying topological constraints. This indicates the decision variable that minimizes the objective function within the curly braces. The value; the value within the curly braces The weighted summation of the deformation energy function and the topological equivalence penalty ensures that the optimization process simultaneously pursues the dual objectives of data alignment accuracy and deformation topology effectiveness; constraints Indicates a grid Each triangular facet in Its Jacobian determinant must be strictly greater than the tolerance threshold. .

[0163] S3.8. Iteratively update the coordinate vector of the second vertex using the Gauss-Newton method. When a certain iteration causes any mesh cell to... Approaching At that time, logarithmic obstacle penalty It tends towards positive infinity, thus using a deterministic gradient backoff mechanism to prevent any parameter updates that might disrupt the homeomorphism.

[0164] Among them, the Gauss-Newton method is a classic nonlinear least squares optimization algorithm. It calculates the optimal update direction and step size for each step by constructing and solving the linearized normal equations through a local quadratic approximation of the objective function in the neighborhood of the current solution.

[0165] A deterministic gradient backoff mechanism refers to a situation where, during optimization iterations, a triangular facet tends to move in the folding direction due to data alignment terms (i.e., its Jacobian determinant begins to approach the tolerance threshold). When, for several terms It will rapidly approach negative infinity, leading to The value instantly soared to positive infinity, correspondingly, The gradient components with respect to the vertex coordinates of the patch will also increase dramatically and point away from the fold boundary. This dramatically increased gradient in the opposite direction is equivalent to erecting an infinitely high energy barrier at the fold boundary—no matter how much driving force the data alignment term generates in trying to push the patch toward the fold, the logarithmic barrier term will generate a greater reverse repulsive force to push the patch away from the fold boundary.

[0166] Figure 6 This diagram illustrates the barrier effect curve in this embodiment, where the logarithmic barrier penalty energy approaches the tolerance threshold as the Jacobian determinant nears its limit. The X-axis represents the Jacobian determinant, and the Y-axis represents the logarithmic barrier penalty energy. Key point: When... Much larger The energy of the penalty term is flat and close to zero when Approaching The energy of the penalty term surges rapidly to positive infinity, forming an energy barrier.

[0167] It is understandable that in the above global optimization model, there is a logarithmic obstacle penalty term. With hard constraints This constitutes a dual topology protection mechanism. This underlying deterministic gradient backoff mechanism forms an absolutely rigid mathematical barrier, fundamentally and completely eliminating the large distortion artifacts such as mesh folding and texture wrinkling that frequently occur under large parallax conditions in traditional registration techniques. It ensures that no matter how large the deformation amplitude, every local region in the registration result maintains a geometrically reasonable and topologically effective homeomorphic mapping relationship.

[0168] Step 4: Calculate the overlapping area between the aligned first and second images, generate image stitching seams based on the dynamic topology seam algorithm, and output the stitched and fused target image.

[0169] The overlapping region refers to the pixel area that is spatially covered by the aligned first and second images, that is, the intersection area where both images have valid pixel values.

[0170] The Dynamic Topology Seam Algorithm is a structure-aware seam planning method that combines the traditional pixel color difference cost with the topological skeleton repulsion field. When planning the seam path, the algorithm not only considers minimizing pixel differences, but also actively avoids the main structural areas in the image.

[0171] An image stitching seam refers to an optimal path determined in the overlapping area, which serves as the boundary line between the pixel sources of the two images in the final stitched result.

[0172] The target image refers to the final panoramic stitching and fusion result synthesized by taking the pixels on both sides of the stitching seam from the aligned first and second images respectively.

[0173] In this embodiment, the most common method in traditional image stitching seam planning techniques is to use the color difference between pixels as the cost indicator and use a dynamic programming algorithm to find the path with the minimum cumulative color difference cost as the stitching seam. However, when facing real-world scenarios with complex subject structures (such as pedestrians, buildings, trees, and other independent objects in the stitching area), the pure color difference method has a serious inherent flaw: when a subject object happens to be located in the overlapping area of ​​two images, and the grayscale distribution of the object in the two images is relatively similar, the dynamic programming algorithm driven by pure color difference is likely to choose a path that directly passes through the middle of the object as the optimal seam, cutting a complete object in half from the middle, forming a clear seam running through the object in the final stitching result, producing an extremely abrupt visual discontinuity effect.

[0174] In this embodiment, the step of generating image stitching seams based on the dynamic topology seam algorithm in S4 includes:

[0175] S4.1 Calculate the pixel-level difference between the first and second images within the overlapping region and generate a difference energy map.

[0176] The pixel-level difference value refers to the Euclidean norm of the difference between the pixel values ​​of the aligned first and second images at each pixel location in the overlapping region.

[0177] A difference energy map is a two-dimensional scalar field formed by arranging the pixel-level difference values ​​at each pixel location within an overlapping region according to their spatial positions. The value at each location represents the degree of visual difference between the two images at that point.

[0178] S4.2. Introduce the pure topological skeleton coordinate system extracted in step 2, and calculate the repulsion weight of each pixel in the overlapping area to the nearest topological skeleton node. The repulsion weight is inversely proportional to the Euclidean distance to the skeleton node.

[0179] Here, the pure topological skeleton coordinate system refers to the pure topological skeleton obtained after trimming in step 2. Spatial coordinate information.

[0180] Repulsion weight is an amplification factor calculated based on the distance from the current pixel position to the nearest topological skeleton point. This weight is inversely proportional to the Euclidean distance to the skeleton node; that is, the closer a pixel position is to the main structural skeleton, the greater its repulsion weight, and the more drastically the cost is amplified. However, the influence of the repulsion field in the background gap region far away from the main structure can be ignored.

[0181] It is understandable that by introducing a repulsive force weight that is inversely proportional to the distance from the topological skeleton, an energy highland or repulsive force field is formed around the main structural skeleton, thereby forcing the splicing seam to actively bypass the main structural region in the subsequent dynamic programming solution.

[0182] S4.3. Construct a cost matrix by multiplying the difference energy map with the repulsion weights. Use a dynamic programming algorithm to find the minimum energy path from the top to the bottom of the image in the cost matrix as the image stitching seam. Force the stitching seam to bypass the topological skeleton region representing the main structure to avoid visual discontinuity.

[0183] The cost matrix refers to the joint cost field that integrates pixel color difference information and topological skeleton repulsion information.

[0184] Dynamic programming is a global path search method based on the Bellman optimality principle. It efficiently solves the globally optimal path that runs through the upper and lower boundaries of an overlapping region through bottom-up recursive calculation.

[0185] The minimum energy path is the continuous path in the cost matrix that minimizes the sum of the joint cost values ​​of all pixel positions along the path.

[0186] Specifically, let the overlapping area formed by the two images after geometric registration in step S3 be . .make Coordinates of the first image after mesh deformation mapping in the overlapping region Pixel value at that location, This represents the pixel value at the same coordinate position in the second image. For each pixel position in the overlapping region... First, the pixel color difference norm of the two images at that location is calculated. Then, the color difference norm is multiplied by a structural repulsion amplification factor controlled by the topological skeleton distance, thereby obtaining the joint cost value that integrates color difference information and structural avoidance information.

[0187] For example, in this embodiment, the joint cost matrix of the fused repulsive force field can be calculated by the following formula:

[0188]

[0189] in, Overlapping region Mid-pixel coordinates The joint cost value at that location is such that a higher value means a greater penalty for placing a seam at that location, and the dynamic programming algorithm will try to avoid this location. For the two images in position The pixel color difference norm at that location, i.e., the value at that position in the difference energy map; (The part in parentheses) This is the structural repulsion amplification factor, and its value is always greater than or equal to 1. This is the topological repulsion coefficient, which controls the overall strength of the repulsive force field of the structure. In practical applications, this coefficient can be adjusted according to the density and importance of the main objects in the scene. For pixels To the set of pure topological skeletons The nearest skeleton point The square of the Euclidean distance, which is located in the denominator, means that when the pixel The closer to the topological skeleton, the larger the value of the fraction and the higher the amplification factor. It is a preset positive smoothing constant, which prevents the denominator from overflowing when the pixel is exactly on the topological skeleton, and at the same time limits the peak value of the repulsive force field to a limited range directly above the skeleton, thereby ensuring the stability of numerical calculation.

[0190] Figure 7 This embodiment illustrates different topological repulsion coefficients. The diagram below shows the distribution of the distance between the seam path and the skeleton. Figure 7 The X-axis represents the topological repulsion coefficient. The value of is given by the Y-axis, which represents the average distance from all pixels along the optimal seam path to the nearest skeleton point. The larger the size, the farther the seam is pushed away from the frame, and the larger the safety margin; at the same time, displaying excessively large... This will cause the seams to be forced to take an excessively long detour, increasing the overall cost of color difference.

[0191] After constructing the aforementioned joint cost matrix Then, the globally optimal seam path that runs through the upper and lower boundaries of the overlapping region can be solved using the classic dynamic programming method.

[0192] For example, in this embodiment, the globally optimal seam path It can be defined by the following formula:

[0193]

[0194] in, The globally optimal seam path is the path that minimizes the sum of the joint costs along the path among all possible continuous paths that traverse the upper and lower boundaries of the overlapping area. This path will serve as the boundary line between the pixel sources of the two images in the final stitching result. This is a variable representing the candidate paths searched during the traversal. To follow the candidate path The sum of the joint values ​​of all pixel locations traversed.

[0195] The solution to the above-mentioned globally optimal path can be efficiently achieved using a bottom-up Bellman dynamic programming recurrence relation, avoiding a brute-force exhaustive search of an exponentially large number of candidate paths.

[0196] For example, in this embodiment, the cumulative energy recursive equation for dynamic programming can be shown as follows:

[0197]

[0198] in, The cumulative energy function is a dynamic programming expression that represents the energy level from the initial boundary of the overlapping region to the pixel position. The minimum cumulative sum of joint cost values ​​among all possible partial paths. The calculation proceeds recursively from the first row downwards, with the cumulative energy value of the first row directly equal to the joint cost value of that row. ; Current pixel position The combined value of the place; To select the value with the minimum accumulated energy from the three adjacent pixel positions in the previous row, these three adjacent positions correspond to the values ​​starting from the top left. Directly above and the upper right The extended path, taking its minimum value, ensures that each step selects the locally optimal sub-path to the current position, conforming to the Bellman optimality principle. When the recursive calculation reaches the last row of the overlapping region, it finds... By identifying the pixel position with the smallest value, and then tracing back from that position along the recorded optimal preceding direction from bottom to top, the complete globally optimal seam path can be reconstructed. .

[0199] It is understandable that in the aforementioned joint cost matrix In the formula, the denominator of the structural repulsion amplification factor includes the square of the shortest distance from the pixel to the topological skeleton. When a candidate path in dynamic programming attempts to approach or traverse the structural skeleton of a main object in the image, the distance between pixels near the skeleton and the skeleton approaches zero, causing a sharp increase in the value of the repulsive force field fraction, which in turn exponentially amplifies the joint cost at that location. In the background gaps between main objects, the influence of the repulsive force field is negligible due to the distance from all skeletons, and the cost is mainly dominated by the color difference term, remaining at a low level. Thus, in the cost matrix of the entire overlapping region, high energy walls are formed near the main structural skeleton, while low energy valleys are formed in the background gaps. Driven by the pursuit of minimizing the global cumulative cost, the dynamic programming algorithm is blocked by these high energy walls and forced to plan a detour path along the low energy valleys in the background gaps. This mathematically guarantees that the optimal seam path will never penetrate any independent main object, fundamentally eliminating the visual fragmentation problem frequently occurring in traditional pure color difference methods, and ensuring that each main object in the final panoramic stitching result maintains a complete, natural, and uninterrupted visual effect.

[0200] Example 2

[0201] This embodiment discloses a graphics image stitching and fusion system based on topology analysis. Specifically, this system can be integrated into an electronic device, such as a terminal or server. The terminal can be a mobile phone, tablet computer, smart Bluetooth device, laptop computer, or personal computer; the server can be a single server or a server cluster composed of multiple servers.

[0202] In this embodiment, the topology-based graphic image stitching and fusion device can also be integrated into multiple electronic devices. For example, the topology-based graphic image stitching and fusion device can be integrated into multiple servers, and the topology-based graphic image stitching and fusion method of this application can be implemented by multiple servers.

[0203] In this embodiment, the server can also be implemented in the form of a terminal.

[0204] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A graphic image stitching and fusion method based on topological analysis, characterized in that, The graphic image stitching and fusion method includes: A first image and a second image are acquired, and the first image and the second image are binarized and partitioned according to a preset grayscale threshold matrix to generate a first binary image and a second binary image. The first binary image and the second binary image are respectively subjected to topological feature parsing. By extracting the topological skeleton of the foreground region of the image and calculating the topological characteristic number, a deterministic topological feature descriptor is generated. Based on the deterministic topological feature descriptor, an initial set of matching point pairs is generated. Using the initial set of matching point pairs as the data driving term, the pixel space of the first image is divided into a two-dimensional non-overlapping triangular mesh set, and a deformation energy function is constructed. During the optimization process of the deformation energy function, a Jacobian determinant positive definite constraint is injected to control the mapping process of the mesh deformation field to satisfy the homeomorphic mapping condition, and the aligned first image is generated. The overlapping region between the aligned first image and the second image is calculated. A cost matrix is ​​constructed within the overlapping region based on pixel-level difference values ​​and repulsion weights of the topological skeleton. The image stitching seam is solved in the cost matrix using a dynamic programming algorithm, and the stitched and fused target image is output. Generate deterministic topological feature descriptors, including: Extract the branch endpoints, intersection nodes, and distance transformation field radius values ​​corresponding to each node from the initial topology skeleton; The coordinates of the branch endpoints, the coordinates of the intersection nodes, the radius values ​​of the distance transformation field of each node, the 0th order Betti number, and the 1st order Betti number are concatenated and combined to generate the deterministic topological feature descriptor containing structural morphology information and topological invariant information. The preset grayscale threshold matrix is ​​shown in the following formula: ； in, This is a grayscale threshold matrix. For the first Each square area The local mean of the grayscale values ​​of all pixels within the range. For the first Each square area The local variance of all pixel grayscale values ​​within the range This is the preset contrast tolerance constant. This is the Heaviside step function, defined as follows: when the value inside the parentheses is positive, the function output is always equal to 1; when the value inside the parentheses is zero or negative, the function output is always equal to 0. This is the background baseline value.

2. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The binarization partitioning process includes: The original image is divided into multiple square blocks, and the mean and variance of grayscale values ​​for each square block are calculated. When the variance of the square block is greater than the preset contrast tolerance, the binarization threshold of the square block is set to the grayscale mean. When the variance of the square block is not greater than the preset contrast tolerance, the binarization threshold of the square block is set to the preset background baseline value. The grayscale value of each pixel in the original image is compared with the binarization threshold corresponding to its square block. Pixels with grayscale values ​​higher than the binarization threshold are identified as foreground structure pixels, and pixels with grayscale values ​​not higher than the binarization threshold are identified as background pixels.

3. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The step of extracting the topological skeleton of the foreground region of the image includes: A median transformation is performed on the first binary image and the second binary image, and a distance transformation field is generated by calculating the Euclidean distance from each foreground pixel to the nearest background boundary pixel within the image. Connect the foreground pixels in the distance transformation field that have at least two equidistant nearest neighbor boundary points on the background boundary to form an initial topological skeleton.

4. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The calculation of the topological characteristic number includes: The initial topological skeleton is abstracted as a one-dimensional simple complex. The intersection nodes and branch endpoints on the initial topological skeleton are abstracted as vertices, and the skeleton segments connecting adjacent vertices are abstracted as edges. Using the Euler-Poincaré formula, the Euler characteristic number, represented by the difference between the 0th and 1st order Betti numbers, is calculated based on the difference between the number of vertices and the number of edges of the simple complex. Based on the boundary operator of the simple complex, a homology group is constructed. The 0th order Betti number and the 1st order Betti number are solved by calculating the rank of each order homology group, where the 0th order Betti number represents the number of connected components and the 1st order Betti number represents the number of independent closed holes.

5. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The deformation energy function includes a feature alignment data term and a mesh rigidity preservation regularization term; The feature alignment data item is calculated in the following way: For each pair of corresponding points in the initial matching point pair set, the centroid coordinates are interpolated according to the coordinates of the second vertex of the grid cell where the first image feature point is located in the corresponding point to obtain the mapping position of the first image feature point, the square Euclidean distance between the mapping position and the position of the second image feature point in the corresponding point is calculated, and the square Euclidean distances of all corresponding points are summed. The mesh rigidity preservation regularization term is calculated as follows: for each pair of adjacent vertices in the two-dimensional non-overlapping triangular mesh set, the square Euclidean norm of the deviation between the deformed edge vector and the original edge vector after rotation by the local optimal rotation matrix is ​​calculated, the square Euclidean norm of all adjacent vertex pairs is summed, and multiplied by the regularization weight coefficient.

6. The graphic image stitching and fusion method based on topology analysis according to claim 5, characterized in that, The injected Jacobian determinant positive definiteness constraint includes: For each mesh cell in the two-dimensional non-overlapping triangular mesh set, the Jacobian matrix corresponding to the local affine mapping of the mesh cell from before deformation to after deformation is calculated based on the coordinates of the first vertex before deformation and the coordinates of the second vertex after deformation. Calculate the Jacobian determinant of the Jacobian matrix; It is mandatory that the Jacobian determinant of each grid cell in the two-dimensional non-overlapping triangular mesh set is strictly greater than a preset tolerance constant, wherein the tolerance constant is a positive value greater than zero.

7. The graphic image stitching and fusion method based on topology analysis according to claim 6, characterized in that, The Jacobian determinant is calculated in the following way: The two edge vectors formed by the coordinates of the second vertex of the deformed grid cell are arranged in columns to form a first matrix. The two edge vectors formed by the coordinates of the first vertex of the grid cell before deformation are arranged in columns to form a second matrix. The first matrix and the inverse of the second matrix are multiplied to obtain the Jacobian matrix of the grid cell.

8. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The cost matrix is ​​constructed based on pixel-level difference values ​​and topological skeleton repulsion weights, including: Calculate the Euclidean norm of the pixel value difference between the aligned first image and the second image at each pixel position within the overlapping region, and generate a difference energy map; Introducing the spatial coordinate information of the topological skeleton, the square of the Euclidean distance from each pixel in the overlapping region to the nearest topological skeleton point is calculated. The square of the square is summed with a repulsion force weight smoothing constant as the denominator. The numerator is a preset positive topological repulsion coefficient. The quotient of the numerator and the denominator is calculated. The value is added to the quotient to obtain the repulsion force weight. The cost matrix is ​​generated by multiplying the Euclidean norm at each pixel location in the difference energy map with the repulsive force weight at that location. The preset positive topological repulsion coefficient is used to control the repulsion strength of the topological skeleton to the path of the image stitching seam.

9. The graphic image stitching and fusion method based on topology analysis according to claim 1, characterized in that, The step of solving for the image stitching seam in the cost matrix using a dynamic programming algorithm includes: Starting from the initial boundary row of the overlapping region, the cumulative energy function value of that row is initialized to the cost matrix value of that row; For each pixel position in each subsequent row, the cost matrix value at that position is added to the minimum value of the cumulative energy function values ​​of that position and its left and right adjacent positions in the previous row to obtain the cumulative energy function value at that position. When the recursion reaches the last boundary row of the overlapping region, starting from the pixel position with the smallest accumulated energy function value in that row, the image stitching seam is obtained by backtracking from the last boundary row to the starting boundary row along the recorded optimal forward direction.

10. A graphics and image stitching and fusion system based on topology analysis, characterized in that, The graphic image stitching and fusion system includes: processor; The memory stores a computer program that, when executed by a processor, implements the topology-based graphic image stitching and fusion method as described in any one of claims 1 to 9.

Images