Non-rigid image registration method, system and use thereof

Abstract

The present application relates to a non-rigid image registration method and system, in particular to a non-rigid image registration method and system based on triangulation and piecewise affine transformation, and belongs to the field of image processing and analysis. The non-rigid image registration method comprises: splitting a feature image to obtain a triangular image block, performing piecewise affine transformation on the triangular image block, and merging all the triangular image blocks after piecewise affine transformation. The present application combines the triangular image blocks registered in sections, and efficiently and accurately realizes the registration of neural circuit imaging across multiple time courses.

Description

TECHNICAL FIELD

[0001] The present application relates to a non-rigid image registration method and system and use thereof, in particular, a non-rigid image registration method and system and use thereof based on triangulation and piecewise affine transformation, and belongs to the field of image processing and analysis. BACKGROUND

[0002] Exploring the characteristics of neural circuits is the key to understanding the cognitive ability of the brain. Neural circuit imaging technology can monitor the activity of hundreds to thousands of neurons at single neuron resolution over multiple days. Long-term recording of multi-time neural imaging is an important technical means to study the functional mechanism of neural circuits, and accurate and effective non-rigid image registration is crucial for long-term recording of multi-time neural circuit imaging. However, during imaging at different times, due to the positional error of microscope repositioning during experimental operation and the non-rigid deformation of the brain, the fields of view at different times usually cannot be accurately matched.

[0003] Before analyzing multi-time neural circuit imaging data, the same features at different times need to be positionally registered. Image registration is widely used in many image analysis tasks, such as computer vision and medical images. A typical image registration process includes four key steps: feature detection, feature matching, mapping function design, and image transformation. In neural circuit imaging registration, the key image features that play an important role include neurons, axons, dendrites, and blood vessels. Most of the established neural image registration tools are developed to solve the registration problem of small deformation in imaging results at the same time, and the registration of neural circuit imaging across multiple times is still a challenge. Traditional image rigid affine transformation has a perfect theoretical basis, but is not suitable for directly solving the non-rigid deformation problem of the brain. Moreover, multi-day multi-time neural circuit imaging has a large complex deformation, and traditional rigid deformation registration methods and methods suitable for small-range non-rigid deformation registration are not suitable for multi-day multi-time neural circuit imaging registration, and are not suitable for directly solving the non-rigid deformation problem of the brain. SUMMARY

[0004] To improve the above technical problems, the present application provides a non-rigid image registration method and system and use thereof, in particular, a non-rigid image registration method and system and use thereof based on triangulation and piecewise affine transformation.

[0005] According to an embodiment of the present application, the non-rigid image registration method comprises: dividing a feature image to obtain triangular image blocks, performing piecewise affine transformation on the triangular image blocks, and merging all the triangular image blocks after piecewise affine transformation.

[0006] According to an embodiment of the present application, the non-rigid image registration method comprises the following steps:

[0007] S1: extracting a feature image from a multi-time neural circuit imaging result;

[0008] S2: identifying a local feature point from the feature image;

[0009] S3: pairing the local feature points;

[0010] S6: expanding the image outwardly to increase the boundary feature points of the expanded image;

[0011] S7: performing a triangular subdivision on the feature image according to the paired feature points to obtain a triangular image block;

[0012] S8: performing a piecewise affine transformation on each triangular image block;

[0013] S9: merging all the triangular image blocks after the piecewise affine transformation to obtain a final non-rigid image pairing result.

[0014] According to an embodiment of the present application, steps S2 and / or S3 of the registration method can be performed automatically or non-automatically.

[0015] According to an embodiment of the present application, when step S2 of the registration method is performed automatically, the registration method can further comprise the following step S4:

[0016] S4: manually adjusting the local feature points, including deleting or adding the local feature points.

[0017] According to an embodiment of the present application, when step S3 of the registration method is performed automatically, the registration method can further comprise the following step S5:

[0018] S5: manually adjusting the local feature point pairs, including deleting or adding the local feature point pairs.

[0019] According to an embodiment of the present application, step S1 can be applicable to a plurality of different neural circuit imaging modes. For example, the imaging mode of the neural circuit can be selected from a desktop two-photon imaging, a small-sized head-mounted two-photon imaging, a small-sized head-mounted single-photon imaging, a three-photon imaging, etc.

[0020] Preferably, step S1 further comprises converting the neuron activity into an optical signal. The way of converting the neuron activity into the optical signal can be selected from OGB, Ga520, virus, transgene, etc.

[0021] Preferably, the object targeted by step S1 can be an animal body, such as selected from mouse, rat, monkey, tree shrew, etc. As an example, the same animal body can be observed multiple times in a multi-time course manner, including multiple observations in different time periods within a day, and / or observations in multiple days. Among them, the imaging results of the neural circuit in each time course can be obtained as an image sequence, and a feature image can be extracted from the imaging image sequence of each time course.

[0022] According to the embodiments of the present application, the way of extracting the feature image from the image sequence can be selected from one or more of the following: Z-axis average, Z-axis maximum projection, Z-axis variance, etc.

[0023] According to the embodiments of the present application, for the feature images of two time courses, one of the feature images can be taken as a static image IM, and the other image can be taken as a dynamic image im. In the image registration process, the static image remains unchanged, and the dynamic image is non-rigidly transformed to the static image to achieve image registration.

[0024] According to the embodiments of the present application, in the process of identifying local feature points from the feature image in step S2, the key local image features include but are not limited to local image features selected from neurons, axons, dendrites, and blood vessels, etc.

[0025] For example, the local feature points can be identified from the feature image in multiple ways, including one or more selected from the following: BRISK feature, MSER feature, ORB feature, SURF feature, KAZE feature, cornerPoints feature based on FAST algorithm, cornerPoints feature based on Harris-Stephens algorithm, cornerPoints feature based on minimum eigenvalue algorithm, etc.

[0026] According to the embodiments of the present application, the better local image features should have rotation invariance, scale invariance, and robustness to noise, etc.

[0027] According to the embodiments of the present application, in the process of local feature point pairing in step S3, especially in the process of automatic pairing, the feature point pairing can be achieved by using multi-scale two-dimensional Fourier transform (2D FFT). Among them, the feature point in the static image is represented as P k =(X k ,Y k ), and the corresponding feature point in the dynamic image is represented as p k . Preferably, P k and p k are in different images, but represent the same object, because the image deformation can cause the spatial positions of the two to be different.

[0028] According to the embodiments of the present application, step S3 can include steps S3-1, S3-2, S3-3 and S3-4:

[0029] S3-1: setting primary two-dimensional Fourier transform feature point pairing parameters;

[0030] S3-2: extracting image block pairs, calculating cross-correlation matrix, and screening potential feature point pairs;

[0031] S3-3: setting secondary two-dimensional Fourier transform feature point pairing parameters;

[0032] S3-4: extracting image block pairs, calculating cross-correlation matrix, and determining final feature point pairs.

[0033] According to the embodiments of the present application, in the process of setting primary two-dimensional Fourier transform feature point pairing parameters in step S3-1, the value S of image block size and the potential pairing feature point cross-correlation threshold CorrThre1 can be set as needed.

[0034] According to the embodiments of the present application, in the process of extracting image block pairs, calculating cross-correlation matrix, and screening potential feature point pairs in step S3-2, the image blocks IM k and im k are extracted from the static and dynamic images respectively with P k1 and p k1 as the center points, and the size of the image blocks is S×S:

[0035] IM k1 (x,y)=IM(X k +x-S / 2,Y k +y-S / 2)

[0036] im k1 (x,y)=im(X k +x-S / 2,Y k +y-S / 2)

[0037] Wherein x,y=1,2,…,S;

[0038] The cross-correlation matrix R1 between the image blocks IM k1 and im k1 is calculated, and the difference DiffTMax1 between the maximum peak value and the second largest peak value of the cross-correlation matrix R1 is calculated;

[0039] The relationship between the difference DiffTMax1 and the pairing feature point cross-correlation threshold CorrThre1 is compared:

[0040]

[0041] If Flag1 = 0, the feature point P k is judged as not paired; k

[0042] If Flag1 = 1, the feature point P k is judged as a potential paired feature point, and then the offset (Δx1, Δy1) of the maximum peak distance 0 point of the cross-correlation matrix R1 is calculated, and the position of the feature point p k in the dynamic image is adjusted according to the offset (Δx1, Δy1); k

[0043] p k = (X k - Δx1, Y k - Δy1).

[0044] According to the embodiment of the present application, in the process of setting the secondary two-dimensional Fourier transform feature point pairing parameters in step S3-3, the value T of the image block size and the potential paired feature point cross-correlation threshold CorrThre2 need to be set;

[0045] As an example, the image block size S in step S3-1 is greater than the image block size T in step S3-3, and as a preferred, the value range of T / S is [0.4 0.8].

[0046] According to the embodiment of the present application, in the process of extracting the image block pair, calculating the cross-correlation matrix, and determining the final feature point pair in step S3-4, the fine-scale 2D FFT is carried out.

[0047] The image blocks IM k and im k are extracted from the static and dynamic images respectively with P k2 and the updated p k2 as the center point, and the size of the image block is T x T:

[0048] IM k2 (x, y) = IM(X k + x - T / 2, Y k + y - T / 2)

[0049] im k2 (x, y) = im(X k + x - Δx1 - T / 2, Y k + y - Δy1 - T / 2)

[0050] Where x, y = 1, 2, …, T;

[0051] The image blocks IM k2 and im k2 ​​Calculate the difference between the maximum peak and the second largest peak of the cross-correlation matrix R2;

[0052] Compare the relationship between the difference DiffTMax2 and the cross-correlation threshold CorrThre2 of paired feature points:

[0053]

[0054] If Flag2 = 0, feature point P k With p k Feature points that are judged as unpaired;

[0055] If Flag2 = 1, feature point P k With p k Once a pair of feature points are identified, the offset (Δx2, Δy2) from the maximum peak of the cross-correlation matrix R2 to point 0 is calculated. Based on the offsets (Δx1, Δy1) and (Δx2, Δy2), the feature point p in the dynamic image is determined. k The final position is denoted as

[0056]

[0057] x k =X k -Δx1-Δx2

[0058] y k =Y k -Δy1-Δy2

[0059] Feature points P located in the static and dynamic images respectively k and This represents the final paired feature point pair.

[0060] According to an embodiment of the present invention, in the process of manually adjusting local feature points in step S4, including deleting and adding local feature points, local feature points can be manually added and deleted through human-computer interaction.

[0061] According to an embodiment of the present invention, in the process of manually adjusting local feature point pairs in step S5, including deleting and adding local feature point pairs, the added local feature points can be manually paired and existing local feature point pairs can be deleted through human-computer interaction.

[0062] According to the embodiment of the present application, in the process of expanding the image outward in step S6 and increasing the feature points on the boundary of the expanded image, the image is expanded outward with the original static image and dynamic image as the center, and the expansion parameter is z, the size of the original image is M*N, and the size of the expanded image is (M*z) x (N*z); as a preferred, the value of the expansion parameter z is in the range of [1.1, 1.5]. The expanded image part can be filled in various ways, such as all 0 filling, all 1 filling, arbitrary given value filling, etc.

[0063] Preferably, on the boundary of the expanded image, the feature points are set in an equidistant manner.

[0064] Preferably, the feature points on the boundary of the expanded static image and dynamic image have the same spatial position.

[0065] According to the embodiment of the present application, in the process of triangulating the feature image according to the paired feature points in step S7, the static image and the dynamic image are triangulated by using the Delaunay triangulation method.

[0066] Preferably, in the process of triangulating the static image and the dynamic image, they have the same triangular grid corresponding order; that is, for any one triangular image block after triangulation in the static image, there is another triangular image block in the triangulation result of the dynamic image, and the two triangular image blocks correspond to three pairs of paired feature points; the static image IM and the dynamic image im are decomposed into N triangular image blocks:

[0067]

[0068]

[0069] wherein imPatch n and imPatch n represent the nth triangular image block in the static and dynamic images, respectively.

[0070] According to the embodiment of the present application, in the process of segmenting and affine transforming each triangular image block in step S8, all the IMPatch n and imPatch n triangular image block pairs are registered. The specific operation is that imPatch n is registered to IMPatch n , and the registered imPatch n is recorded as imPatchReg n .

[0071] imPatch nto IMPatch n The transformation matrix of registration is denoted as TF:

[0072]

[0073] For each pair of feature points, P k = (X k , Y k ),

[0074] [X k Y k 1] = [x k y k 1] · TF

[0075]

[0076] Preferably, each group of triangular image patches imPatch n to IMPatch n has 3 pairs of feature points, where the feature points of imPatch n are denoted as (x 1n , y 1n ), (x 2n , y 2n ), (x 3n , y 3n ), and the feature points of IMPatch n are denoted as (X 1n , Y 1n ), (X 2n , Y 2n ), (X 3n , Y 3n ):

[0077]

[0078] Preferably, the transformation matrix TF is solved by solving the above equation set, and the registered image patch imPatchReg n is obtained according to the affine transformation.

[0079] imPatchReg n = warp(imPatch n , TF)

[0080] According to the embodiment of the present application, in step S9, all the triangular image patches after segmented affine transformation are merged, and in the process of non-rigid image registration, all the registered triangular image patches imPatchReg n are combined to obtain the registered dynamic image imReg:

[0081]

[0082] According to the embodiments of the present application, there is also provided a non-rigid image registration system for implementing the registration method, comprising an image input device, an image extraction and feature point identification device, an image partitioning device, a partitioned image processing device and a result output device suitable for the registration method.

[0083] The present application also provides the use of the registration method or the system for cross-time course neural circuit imaging registration.

[0084] Advantages

[0085] The technical solution of the present application is a non-rigid image registration method and system based on triangular partitioning and piecewise affine transformation, which realizes feature point pairing in a multi-scale two-dimensional Fourier transform (2D FFT) manner, performs triangular partitioning on static images and dynamic images according to the paired feature points, and carries out piecewise affine transformation on the triangular image blocks after triangular partitioning, thereby realizing image registration. Moreover, the triangular image blocks after piecewise registration are combined, thereby efficiently and accurately realizing cross-time course neural circuit imaging registration. BRIEF DESCRIPTION OF DRAWINGS

[0086] Figure 1 A schematic diagram of the steps of the non-rigid image registration method based on triangular partitioning and piecewise affine transformation of the present application;

[0087] Figure 2 A schematic diagram of steps 3-1 to 3-4 of the non-rigid image registration method of the present application;

[0088] Figure 3 Results of image registration on the multi-time course neural circuit imaging results with spatial position deformation in Example 1;

[0089] Figure 4 Results of local feature point identification from static images and dynamic images in Example 2;

[0090] Figure 5 Feature point pairing identification effect based on two-dimensional Fourier transform in Example 3;

[0091] Figure 6 Results of piecewise affine transformation on the triangular image blocks after triangular partitioning and merging of all the triangular image blocks after piecewise affine transformation to complete non-rigid image registration in Example 4;

[0092] Figure 7The result of the non-rigid image registration of the triangular image blocks after the segmentation affine transformation and the merging of all the triangular image blocks after the segmentation affine transformation in Example 5 is shown in FIG. 6.

[0093] Figure 8 The effect of the image registration before and after in Example 6 is shown in FIG. 7. DETAILED DESCRIPTION

[0094] The technical solutions of the present application will be further described in detail below in combination with specific examples. It should be understood that the following examples are only illustratively described and explained, and should not be interpreted as limiting the scope of protection of the present application. Any technology achieved based on the above description of the present application is covered within the scope of protection intended by the present application.

[0095] Unless otherwise specified, the instruments used in the following examples are commercially available.

[0096] Example 1

[0097] This example provides a non-rigid image registration method and system based on triangulation and segmentation affine transformation, which is used for image registration of multi-time neural circuit imaging results with spatial position deformation.

[0098] In this example, two-photon (2p) neural circuit imaging technology was used to record the neural circuit activity of the PFC cortex of mice. All animal protocols were approved by the Institutional Animal Care and Use Committee (IACUC) of Peking University. The experimental procedures were carried out in accordance with the Guide for the Care and Use of Laboratory Animals. All mice were housed in the animal facility of the Peking University Experimental Animal Center at 25±2℃, with a 12-hour light / dark cycle until used for 2p imaging experiments. A two-photon microscope (B-Scope, Thor Labs) and a femtosecond laser (MaiTai BB DS-OL, Spectra-Physics) were used to monitor the fluorescence transients of cortical neurons in the prefrontal cortex (mPFC) of mice. The laser wavelength was set to 920 nm. A 16x water immersion objective lens (0.8NA, Nikon) was used, with an image field of 517.77x517.77μm and a pixel of 768x768. The scanning frequency was 20.4Hz. AAV9.syn.GCaMP6f.WPRE.SV40 (Penn Vector Core) was used to label the calcium signals of cortical neurons. The 2p calcium imaging experiment started 7 days after the head fixation surgery, which was recorded as day 0; the time courses 1-5 of 2p imaging recording were performed on days 0, 3, 6, 9 and 12. Each time course recorded 1300 frames of image sequences. Penn Vector Core) was used to label the calcium signals of cortical neurons. The 2p calcium imaging experiment started 7 days after the head fixation surgery, which was recorded as day 0; the time courses 1-5 of 2p imaging recording were performed on days 0, 3, 6, 9 and 12. Each time course recorded 1300 frames of image sequences.

[0099] As Figure 3As shown in the figure, among them, the static image (Session 1) is the characteristic image of the first time course neural circuit image sequence, the dynamic image (Session N) is the characteristic image of the Nth time course neural circuit image sequence, and the fusion image is the result of fusing the static image and the dynamic image together. As can be seen from the fusion image, there is a large misalignment between the reference image and the target image.

[0100] Embodiment 2

[0101] This embodiment provides a non-rigid image registration method and system based on triangulation and piecewise affine transformation according to the present application, which identifies the results of local feature points from static images and dynamic images.

[0102] As Figure 4 shown in the figure, in this embodiment, the feature point recognition method adopts the MSER (Maximally Stable Extremal Regions) feature recognition. The basic principle of MSER is to take the threshold value of a gray image (gray value is 0-255) for binaryzation processing, and the threshold value increases from 0 to 255. The increase of the threshold value is similar to the rise of the water surface in the watershed algorithm. With the rise of the water surface, some lower hills will be submerged, and if you look down from the sky, the earth will be divided into land and water areas, which is similar to a binary image. In all the binary images obtained, some connected regions in the image change very little, or even not at all, and then this region is called the maximum stable extremal region. This is similar to when the water surface continues to rise, the area of some places submerged by water does not change. MSER has the following characteristics: affine transformation invariance to image gray scale; stability, only the regions supported within the same threshold range will be selected; multi-scale detection can be achieved without any smoothing processing, i.e. small and large structures can be detected.

[0103] In this embodiment, MSER feature recognition is realized based on the MATLAB R2021a software platform, and the detectMSERFeatures function is used to realize feature recognition, with the parameters set as RegionAreaRange: [40 100]; MaxAreaVariation: 0.1; ThresholdDelta: 1. The detected original MSER features are further screened to obtain the final MSER features; the screening method is to calculate the Intersection over Union (IoU) between all detected MSER features, and for two MSER features with an Intersection over Union greater than a threshold value (0.2), the MSER feature with a larger area is retained and the MSER feature with a smaller area is deleted.

[0104] Embodiment 3

[0105] This embodiment provides a non-rigid image registration method and system based on triangulation and piecewise affine transformation as described in this invention, which achieves feature point pairing and recognition based on two-dimensional Fourier transform.

[0106] like Figure 5 As shown, in this embodiment, the pairing parameters for the first-order two-dimensional Fourier transform feature points are: the image patch size S = 36; and the potential cross-correlation threshold CorrThre1 = 0.5. For the three MSER feature points shown in the figure, image patches P1, P2, and P3 are extracted from the static image, and their corresponding image patches p1, p2, and p3 in the dynamic image are extracted. The two-dimensional fast Fourier transform is implemented using the fft2 function on the MATLAB R2021a software platform. The cross-correlation matrix between the image patch pairs is calculated, and the difference between the maximum and second-largest peak values ​​in the cross-correlation matrix, DiffTMax, is calculated. For images I1 and I2, the cross-correlation matrix is ​​calculated as follows:

[0107] FI1 = fft2(I1);

[0108] FI2 = fft2(I2);

[0109] FR=FI1.*conj(FI2);%calculating correlation

[0110] R = ifft2(FR);

[0111] R = fftshift(R);

[0112] The DiffTMax values ​​for image patch pairs P1 and p1, P2 and p2, and P3 and p3 are 3.45, 1.30, and 0.39, respectively. Based on the set cross-correlation threshold CorrThre1 (0.5), P1 and p1, P2 and p2 are determined to be potential feature point pairs, while P3 and p3 are determined to be unpaired feature points.

[0113] Further, the pairing parameters for the second-order two-dimensional Fourier transform feature points were set as follows: the image patch size was T = 24; the cross-correlation threshold for potential paired feature points was CorrThre2 = 0.5. For potential feature point pairs P1 and p1, P2 and p2, fine-grained feature point pairing screening was performed.

[0114] Example 4

[0115] This embodiment provides a non-rigid image registration method and system based on triangulation and piecewise affine transformation as described in this invention, which results in the triangulation of an image based on paired feature points.

[0116] As Figure 6 shown in the embodiment, in which the original paired feature points in the reference image and the target image are marked with circles and plus signs, respectively. There are a total of 8 pairs of original feature points; the original reference image and the target image are expanded with an expansion parameter z = 1.2, and the expanded area is filled with 0; 8 pairs of feature points are established on the boundary of the expanded image, including 4 pairs of corner points and 4 pairs of image boundary midpoints; Delaunay triangulation is realized based on the MATLAB R2021a software platform, and the specific implementation function is delaunay; a total of 22 triangular image blocks are obtained.

[0117] Embodiment 5

[0118] This embodiment provides a non-rigid image registration method and system based on triangulation and segmented affine transformation, which performs segmented affine transformation on the triangular image blocks after triangulation, and merges all the triangular image blocks after segmented affine transformation, to complete the registration result of non-rigid image registration, as Figure 7 shown.

[0119] In this embodiment, affine transformation is performed on each pair of triangular image blocks obtained after triangulation, as Figure 6 shown. The result of affine transformation is: transform each triangular image block in the dynamic image to the position of the corresponding triangular image block in the static image; that is, map the content of the dynamic image to the corresponding position of the static image, thereby realizing the registration of each triangular image block. Affine transformation of each pair of triangular image blocks is realized based on the MATLAB R2021a software platform, and the specific implementation function is estimateGeometricTransform, with 3 pairs of position coordinates of the triangular image blocks as input and affine as the method; the transformation matrix TF is calculated.

[0120] Embodiment 6

[0121] This embodiment provides a non-rigid image registration method and system based on triangulation and segmented affine transformation, and the effects before and after image registration.

[0122] As Figure 8 shown, before registration, there is a large spatial position misalignment between the image contents in the fusion image of the static image and the dynamic image. After registration, the mismatch between the image contents in the static image and the registered fusion image is eliminated.

[0123] The above has exemplarily described the embodiments of the technical scheme of the present application. It should be understood that the protection scope of the present application is not limited to the above embodiments. Any modification, equivalent replacement, improvement, etc. made by those skilled in the art within the spirit and principle of the present application shall be included in the protection scope of the claims of the present application.

Claims

1. A non-rigid image registration method, wherein the method comprises: The feature image is divided into triangular image blocks, and the triangular image blocks are subjected to piecewise affine transformation. All the triangular image blocks after piecewise affine transformation are then merged. The non-rigid image registration method includes the following steps: S1: Extract feature images from multi-temporal neural circuit imaging results; S2: Identify local feature points from the feature image; S3: Pair up local feature points; S6: Expand the image outward and add feature points to the boundary of the expanded image; S7: Based on the paired feature points, the feature image is triangulated to obtain triangular image blocks; S8: Perform a piecewise affine transformation on each triangular image block; S9: Merge all the triangular image blocks after piecewise affine transformation to obtain the final non-rigid image pairing result; Step S1 targets an animal body; Use a multi-time-course approach to observe a specific neural circuit of the same animal body multiple times, including multiple observations at different time periods within a day, and / or observations over multiple days; Obtain an image sequence from the neural circuit imaging results of each time-course, and extract a feature image from the imaging image sequence of each time-course. For feature images from two time periods, one feature image is designated as the static image IM, and the other as the dynamic image im. During image registration, the static image remains unchanged, while the dynamic image undergoes a non-rigid transformation to become the static image, thus achieving image registration.

2. The registration method as described in claim 1, wherein steps S2 and / or S3 of the registration method are performed automatically or not automatically; When step S2 of the registration method is performed automatically, the registration method further includes the following step S4: S4: Manually adjust local feature points, including deleting and adding local feature points; and / or When step S3 of the registration method is performed automatically, the registration method further includes the following step S5: S5: Manually adjust local feature point pairs, including deleting and adding local feature point pairs.

3. The registration method as described in claim 1, wherein step S1 is applicable to neural circuit imaging modes selected from the following: desktop two-photon imaging, small head-mounted two-photon imaging, small head-mounted single-photon imaging, and three-photon imaging.

4. The registration method as described in claim 1, step S1 further includes converting neuronal activity into light signals, wherein the method of converting neuronal activity into light signals is selected from OGB, viruses, and genetically modified organisms.

5. The registration method as described in claim 1, wherein the method for extracting feature images from the image sequence is selected from one or more of the following: Z-axis averaging, Z-axis maximum projection, and Z-axis variance.

6. The registration method according to any one of claims 1-5, wherein in the process of identifying local feature points from the feature image in step S2, the local image features include, but are not limited to, local image features selected from neurons, axons, dendrites and blood vessels; Local feature points identified from feature images are selected from one or more of the following: BRISK features, MSER features, ORB features, SURF features, KAZE features, cornerPoints features based on the FAST algorithm, cornerPoints features based on the Harris–Stephens algorithm, and cornerPoints features based on the minimum eigenvalue algorithm.

7. The registration method according to any one of claims 1-5, wherein in the local feature point pairing process of step S3, feature point pairing is achieved by using a multi-scale two-dimensional Fourier transform (2D FFT); wherein, Feature points in a static image are represented as P k =(X k Y k The potential corresponding feature points in the dynamic image are represented as p. k .

8. The registration method as described in claim 7, P k With p k They are in different images, but represent the same object.

9. The registration method as described in claim 7, wherein, Step S3 includes the following steps S3-1, S3-2, S3-3, and S3-4: S3-1: Set the pairing parameters for the first-order two-dimensional Fourier transform feature points; S3-2: Extract image patch pairs, calculate cross-correlation matrix, and screen potential feature point pairs; S3-3: Set the pairing parameters for the feature points of the second-order two-dimensional Fourier transform; S3-4: Extract image patch pairs, calculate the cross-correlation matrix, and determine the final feature point pairs.

10. In the registration method as described in claim 9, during the process of setting the pairing parameters of the first-order two-dimensional Fourier transform feature points in step S3-1, the value of the image block size S and the potential cross-correlation threshold CorrThre1 of the paired feature points can be set as needed.

11. The registration method as described in claim 9, in the process of extracting image patch pairs, calculating the cross-correlation matrix, and screening potential feature point pairs in step S3-2, respectively using P k With p k Using the center point, extract image patches (IM) from static and dynamic images. k1 and im k1 The size of the image patch is S×S: IN THE k1 (x,y)<IM(X k +xS / 2,Y k +yS / 2) in the k1 (x,y)<im(X k +xS / 2,Y k +yS / 2) Where x, y = 1, 2, ..., S; X k Y k The center point P is respectively K coordinates (X) k Y k The x-coordinate and y-coordinate values ​​of the point; Calculate image patch IM k1 and im k1 Given the cross-correlation matrix R1, calculate the difference DiffTMax1 between the maximum peak and the second largest peak of the cross-correlation matrix R1. Compare the relationship between the difference DiffTMax1 and the cross-correlation threshold CorrThre1 of paired feature points: If Flag1 = 0, feature point P k With p k Feature points that are judged as unpaired; If Flag1 = 1, feature point P k With p k The feature points are identified as potential paired feature points. Then, the offset (Δx1, Δy1) from the maximum peak of the cross-correlation matrix R1 to the zero point is calculated. Based on the offset (Δx1, Δy1), the feature points p in the dynamic image are adjusted. k Location; p k =(X k -Δx1,Y k -Δy1)。 12. The registration method as described in claim 9, wherein, In step S3-3, when setting the pairing parameters of the second-order two-dimensional Fourier transform feature points, it is necessary to set the value of the image patch size T and the potential cross-correlation threshold CorrThre2 of the paired feature points.

13. The registration method as described in claim 12, wherein the image block size S in step S3-1 is greater than the image block size T in step S3-3, and the value range of T / S is [0.4 0.8].

14. In the registration method as described in claim 11, during the process of extracting image patch pairs, calculating the cross-correlation matrix, and determining the final feature point pairs in steps S3-4, a fine-scale 2DFFT is performed; P respectively k With the updated p k Using the center point, extract image patches (IM) from static and dynamic images. k2 and im k2 The size of the image patch is T×T: IM k2 (x,y)=IM(X k +x-T / 2,Y k +y-T / 2) im k2 (x,y)=im(X k +x-Δx1-T / 2,Y k +y-Δy1-T / 2) Where x, y = 1, 2, ..., T; Calculate image patch IM k2 and im k2 Calculate the difference between the maximum peak and the second largest peak of the cross-correlation matrix R2; Compare the relationship between the difference DiffTMax2 and the cross-correlation threshold CorrThre2 of paired feature points: If Flag2 = 0, feature point P k With p k Feature points that are judged as unpaired; If Flag2 = 1, feature point P k With p k Once a pair of feature points are identified, the offset (Δx2, Δy2) from the maximum peak of the cross-correlation matrix R2 to the zero point is calculated. Based on the offsets (Δx1, Δy1) and (Δx2, Δy2), the feature point p in the dynamic image is determined. k The final position is denoted as x k =X k -Δx1-Δx2 y k =Y k -Δy1-Δy2 Feature points P located in the static and dynamic images respectively k and This is the final paired feature point pair.

15. The registration method according to any one of claims 1-5, wherein in the process of manually adjusting local feature points in step S4, including deleting and adding local feature points, local feature points are manually added and deleted through human-computer interaction.

16. The registration method according to any one of claims 1-5, wherein in step S5, the local feature point pairs are manually adjusted, including the deletion and addition of local feature point pairs, the added local feature points are manually paired and the existing local feature point pairs are deleted through human-computer interaction.

17. The registration method according to any one of claims 1-5, in the process of expanding the image outward and adding feature points of the expanded image boundary in step S6, the image is expanded outward with the original static image and dynamic image as the center, the expansion parameter is z, the size of the original image is M×N, and the size of the expanded image is (M*z)×(N*z).

18. In the registration method as described in claim 17, the range of the extended parameter z is [1, 1, 1, 5].

19. The registration method as described in claim 17, wherein feature points are set at equal intervals on the expanded image boundary; the feature points on the boundaries of the expanded static image and the dynamic image have the same spatial position.

20. The registration method according to any one of claims 1-5, wherein in step S7, during the process of triangulating the feature image according to the paired feature points, the static image and the dynamic image are triangulated using the Delaunay triangulation method.

21. The registration method as described in claim 14, wherein during the triangulation of the static image and the dynamic image, they have the same triangular mesh correspondence order; that is, for any triangular image block after triangulation in the static image, there exists another triangular image block in the result of triangulation in the dynamic image, and these two triangular image blocks correspond to three sets of paired feature point pairs; the static image IM and the dynamic image im are decomposed into N triangular image blocks: IMPatch n and imPatch n These represent the nth triangular image patch in the static and dynamic images, respectively.

22. The registration method as described in claim 21, wherein during step S8, when performing piecewise affine transformation on each triangular image block, all IMPatch... n and imPatch n Registration of triangular image patch pairs: imPatch n To IMPatch n Registration, the registered imPatch n Recorded as imPatchReg n ; imPatch n To IMPatch n The transformation matrix for registration is denoted as TF: For each paired feature point pair, P k =(Xk,Yk), [X k AND k 1]=[x k and k 1]·TF 23. The registration method as described in claim 22, wherein each set of triangular image patches corresponds to imPatch n To IMPatch n It has 3 sets of paired feature points, where imPatch n The feature points are denoted as (x 1n ,y 1n ),(x 2n ,y 2n ),(x 3n ,y 3n IMPatch n The feature points are denoted as (X) 1n ,Y 1n ),(X 2n ,Y 2n ),(X 3n ,Y 3n ): Solving the system of equations yields the transformation matrix TF, and the registered image patch imPatchReg is obtained through affine transformation. n ; imPatchReg n =warp(imPatch n ,TF)。 24. In the registration method as described in claim 22, during step S9, when merging all the segmented affine transformation triangular image blocks to complete the non-rigid image registration, all registered triangular image blocks imPatchReg are... n Combined, we obtain the registered dynamic image imReg:

Images