A CUDA acceleration-based large three-dimensional image morphing method

By segmenting ultra-large 3D images and utilizing a parallel architecture accelerated by GPUs, the problem of image distortion under computer memory constraints was solved, achieving efficient image registration, significantly improving operating efficiency and shortening processing time.

CN115984098BActive Publication Date: 2026-04-24ANHUI UNIV
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
ANHUI UNIV
Filing Date
2023-01-06
Publication Date
2026-04-24

AI Technical Summary

Technical Problem

Existing technologies suffer from excessively long computation times due to computer memory limitations when processing ultra-large 3D images, making it difficult to achieve rapid deformation operations. This is especially true in high-throughput image analysis tasks, where existing tools cannot effectively improve the efficiency of image registration.

Method used

A CUDA-accelerated method is used to divide ultra-large 3D images into blocks, and the parallel architecture of the GPU is used to perform linear interpolation calculations. The affine transformation matrix is ​​used to deform and stitch the image blocks, achieving rapid image deformation and registration.

Benefits of technology

While keeping computer memory limits in mind, it significantly improves the efficiency of image registration, enabling it to process ultra-large 3D image data that is more than twice the size of computer memory, achieving a speedup of approximately 11 times and greatly shortening the time for deformation registration.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN115984098B_ABST
    Figure CN115984098B_ABST
Patent Text Reader

Abstract

The present application relates to a kind of based on CUDA acceleration's super large three-dimensional image morphing method, comprising: the affine transformation matrix of two sets of registration feature points is calculated by affine transformation;Target image is blocked, and the image morphing of each small block image block obtained after blocking is carried out;Utilize affine transformation matrix, realize three-dimensional image voxel traversal parallelization operation by CUDA;Read the 3D image block after morphing, splice into the 2D continuous sequence image of original image size.The present application utilizes the principle of target image blocking, and utilizes the parallel architecture of GPU to accelerate the calculation process of image linear interpolation, while solving can be for super large three-dimensional image morphing, also effectively improves the running efficiency of image registration;Experiment proves that this method can process more than twice the computer memory super large three-dimensional image data under the platform of CPU as i5-10400F, GPU as RTX-3060, memory as 16G, and about 11 times speedup ratio is obtained, greatly shorten the time of morphing registration.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of image registration and parallelization acceleration technology, and in particular to a method for deformation of ultra-large three-dimensional images based on CUDA acceleration. Background Technology

[0002] Transforming images to standard coordinate space is crucial for many image-computation-related tasks. However, for multi-dimensional and high-resolution images, precise deformation operations consume significant computer memory and are computationally very time-consuming. High-throughput image analysis research, such as brain mapping and registration projects, requires high-performance image transformation tools compatible with common image analysis pipelines. With the continuous advancement of modern microscopic imaging techniques, the size of generated biological image data quickly exceeds the size of computer memory, making it difficult to load all of this biological image data into memory, and computer hardware configurations often do not support this. Therefore, methods for processing ultra-large-scale image data are needed.

[0003] Affine transformation is essentially a simple translation combined with some linear transformations. Its key characteristic is "preserving points," "preserving lines," and "preserving surfaces," meaning that the order of points in the transformed image remains unchanged, and parallel features in the image remain parallel. BigDataViewer, a plugin for the image processing software package Fiji, fills the gap in Fiji's capabilities for browsing ultra-large-scale image data. Users can browse TB-level 5D data locally or remotely. It is primarily based on the custom HDF5 (Hierarchical Data Format), a file format designed for storing and organizing large amounts of data. It supports multiple data types and is specifically designed for flexible and efficient I / O and managing large volumes and complex data. Vaa3D expanded its TB-level data processing capabilities in 2017 and integrated it into the TeraFly plugin. TeraFly adopts a pyramid-level hierarchical data organization method and has built-in fast data format conversion; it maps user clicks on the 2D screen to 3D data, allowing users to move the current view simply by clicking on areas of interest, greatly speeding up image browsing.

[0004] When deforming ultra-large 3D images, given the limitations of computer memory, the image must be processed in blocks. However, precise image deformation operations result in extremely long CPU computation times. Furthermore, for ultra-large 3D images, deformation processing needs to be performed on each individual image block, making the entire process incredibly time-consuming, even on powerful workstations. Therefore, finding a fast deformation method for ultra-large 3D images within the constraints of computer memory has become an urgent technical problem to be solved. Summary of the Invention

[0005] The purpose of this invention is to provide a CUDA-accelerated method for deforming ultra-large 3D images while effectively improving the efficiency of image registration.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a CUDA-accelerated method for deforming ultra-large three-dimensional images, the method comprising the following sequential steps:

[0007] (1) Calculate the affine transformation matrix of the two sets of registration feature points through affine transformation;

[0008] (2) Divide the target image into blocks, and perform image deformation on each small image block obtained after the block division to obtain the deformed 3D image block: using the affine transformation matrix, the three-dimensional image voxel traversal parallel operation is implemented through CUDA, and trilinear interpolation is performed on each image block.

[0009] (3) Perform image stitching operation for each small image block: Read the deformed 3D image blocks and stitch them into a 2D continuous sequence of images of the original image size.

[0010] Step (1) specifically refers to: obtaining two sets of registration feature points, each with N points; normalizing the two sets of registration feature points; calculating the affine transformation matrix according to the affine transformation formula; and finally, inverse normalizing the affine transformation matrix to obtain the final affine transformation matrix as follows:

[0011]

[0012] Among them, α 11 , ɑ 12 ,...,a 34 All are affine parameters.

[0013] Step (2) specifically includes the following steps:

[0014] (2a) Divide the target image into blocks to obtain small image blocks;

[0015] (2b) Perform deformation processing on each small image block to obtain the deformed 3D image block;

[0016] (2c) Read the deformed 3D image block into memory, use CUDA to implement the parallel operation of voxel traversal of the three-dimensional image, use GPU-accelerated trilinear interpolation to process any small image block, i.e., image block A, and finally name and save the file according to the offset of image block A, and save it in 3D image format.

[0017] (2d) Process all small image blocks according to the processing procedure of image block A, and finally generate N. x *N y *N z N small image blocks, after being deformed, are placed in a folder. x N y N z This represents the number of blocks that can be taken along the X, Y, and Z axes.

[0018] Step (3) specifically includes the following steps:

[0019] (3a) First, all the small image patches of the target image form a coordinate system, ranging from [0,0,0 to [N]. x -1,N y -1,N z -1], therefore all small image blocks are divided into N z Layered data, each layer has N layers. x *N y There are several small image blocks, and let any layer k∈[0,N] be an example. z The Z-axis size of the small image patch [-1] is Z max The standard block size is set to [X]. b ,Y b Z b And satisfy X b =Y b =Z b =S per Under the condition of Z max The definition is as follows:

[0020]

[0021] Among them, Z n The Z-axis size of the target image;

[0022] (3b) Set a control variable Z during the splicing process. control Control variable Z control Used to control the amount of memory consumed when reading small image blocks;

[0023] (3c) When Z max ≤Z control At that time, read N of the kth layer x *N y Each small image block is read as Z. max For layer data, let the coordinates of any voxel within any small image block, i.e., image block A, be [x, y, z]. Then the voxel coordinates of this voxel in the target image are [x+X]. offset y+Y offset z+Z offset The offset of the origin of the starting coordinates of image block A in the target image is set to [X]. offset Y offset Z offset The image is stitched together using the X and Y axis offsets of each small image block, ultimately generating a Z-shaped image. max Zhang Image;

[0024] (3d) When Z max >Z control At this time, the processing steps of reading and stitching images K times are required, and their definitions are as follows:

[0025]

[0026] Each time, N is read from the k-th layer. x *N y Each small image block is read as Z. control The data from the layer is processed, and the image stitching process mentioned in step (3c) is performed to generate Z. control Zhang slices;

[0027] (3e) N z The layer data are processed in steps (3c) and (3d) to finally generate Z. n Zhang 2D continuous sequence images.

[0028] Step (2a) specifically includes the following steps:

[0029] (2a1) First, obtain the size of the target image. Since the ultra-large 3D image data is stored as a continuous 2D sequence of images in a folder and arranged in Z-axis order, read any image in the folder and obtain its X and Y axis dimensions, denoted as X. n Y n This refers to the X and Y axis dimensions of the target image, and the number of images read from the folder is denoted as Z. n This refers to the Z-axis size of the target image;

[0030] (2a2) Define a standard block size [X] b ,Y bZ b And satisfy X b =Y b =Z b =S per Let N x N y N z The number of blocks that can be taken along the X, Y, and Z axes is defined as follows:

[0031]

[0032]

[0033]

[0034] (2a3) If the target image boundary is encountered, the size of the small image block is [X]. b′ ,Y b′ Z b′ The specific formula is as follows:

[0035] X b′ =X n -N x *X b

[0036] Y b′ =Y n -N y *Y b

[0037] Z b′ =Z n -N z *Z b

[0038] (2a4) The final number of small image blocks is determined to be N based on the target image size. sum Its definition is as follows:

[0039] N sum =N x *N y *N z .

[0040] In step (2b), the image deformation of each small image block specifically includes the following steps:

[0041] (2b1) After the target image is divided into blocks, a coordinate system is set according to the origin of the target image. Each small image block then obtains a corresponding image block number, i.e., the number range is:

[0042] [0,0,0] to [N] x -1,N y -1,Nz -1]

[0043] For any image patch A, let its corresponding image patch number be [a,b,c], and its image patch size be [X]. A ,Y A Z A If ], then its X-axis dimension is X. b or X b′ Its Y-axis dimension is Y b Or Y b′ Its Z-axis dimension is Z b or Z b′ ;

[0044] (2b2) ​​Since each small image block is processed, the offset of the origin of the starting coordinates of image block A in the target image is set to [X]. offset Y offset Z offset ], which is defined as:

[0045] X offset =a*X b

[0046] Y offset =b*Y b

[0047] Z offset =c*Z b

[0048] The X and Y axis offsets of each image block are recorded and stored in two lists respectively;

[0049] (2b3) Calculate the coordinates of the eight vertices of image block A based on the offset obtained in step (2b2);

[0050] (2b4) Use the affine transformation matrix to convert the coordinates of the eight vertices of image block A into the coordinates of the eight vertices of the original image, and sort the coordinates in the three directions of X, Y and Z. Select the point with the largest and smallest values ​​in these three directions. The difference between the two values ​​is the range of the selected image block in the original image.

[0051] Step (2c) specifically includes the following steps:

[0052] (2c1) First, allocate N to the CUDA kernel function. Z ×N Y N threads, of which N Z N represents the Z-axis size of the 3D image to be registered. Y Set the Y-axis size of the 3D image to be registered, and then fix the blockDim in the kernel function to a 2D block with a size of 32*32;

[0053] (2c2) Each thread has a unique row number and column number, with the row number ranging from 0 to N. Z Its column number ranges from 0 to N. Y The size of the Z and Y axes corresponds to that of the three-dimensional image data;

[0054] (2c3) Then put the CUDA kernel function into a space of size 0 to N. X In the loop, N X The X-axis size is the size of the 3D image to be registered. Finally, each thread processes the work of each voxel in the 3D image to be registered. Assuming that this is the kth loop, the thread with thread number (m, n) processes the voxel corresponding to the X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the 3D voxel. That is, each thread number corresponds to one voxel in the 3D image to be registered.

[0055] As can be seen from the above technical solution, the beneficial effects of the present invention are as follows: First, the present invention utilizes the principle of dividing the target image into blocks and uses the parallel architecture of the GPU to accelerate the calculation process of linear interpolation of the image, which can effectively improve the running efficiency of image registration while deforming ultra-large three-dimensional images; Second, experiments have shown that the present method can process ultra-large three-dimensional image data with more than twice the computer memory on a platform with CPU i5-10400F, GPU RTX-3060 and memory 16G, and achieves an acceleration ratio of about 11 times, which greatly shortens the deformation registration time. Attached Figure Description

[0056] Figure 1 This is a flowchart of the method of the present invention;

[0057] Figure 2 This is a schematic diagram of image block A in this invention. Detailed Implementation

[0058] like Figure 1 As shown, a CUDA-accelerated method for deforming ultra-large 3D images includes the following sequential steps:

[0059] (1) Calculate the affine transformation matrix of the two sets of registration feature points through affine transformation;

[0060] (2) Divide the target image into blocks, and perform image deformation on each small image block obtained after the block division to obtain the deformed 3D image block: using the affine transformation matrix, the three-dimensional image voxel traversal parallel operation is implemented through CUDA, and trilinear interpolation is performed on each image block.

[0061] (3) Perform image stitching operation for each small image block: Read the deformed 3D image blocks and stitch them into a 2D continuous sequence of images of the original image size.

[0062] Step (1) specifically refers to: obtaining two sets of registration feature points, each with N points; normalizing the two sets of registration feature points; calculating the affine transformation matrix according to the affine transformation formula; and finally, inverse normalizing the affine transformation matrix to obtain the final affine transformation matrix as follows:

[0063]

[0064] Among them, α 11 , ɑ 12 ,...,a 34 All are affine parameters.

[0065] Step (2) specifically includes the following steps:

[0066] (2a) Divide the target image into blocks to obtain small image blocks;

[0067] (2b) Perform deformation processing on each small image block to obtain the deformed 3D image block;

[0068] (2c) Read the deformed 3D image block into memory, use CUDA to implement the parallel operation of 3D image voxel traversal, use GPU-accelerated trilinear interpolation to process any small image block, i.e., image block A, and finally name and save the file according to the offset of image block A, and save it in the 3D image format "V3DRAW"; CUDA (Compute Unified Device Architecture) is a general-purpose parallel computing architecture launched by NVIDIA. This architecture enables GPUs to solve complex computing problems and implement the parallel operation of 3D image voxel traversal.

[0069] (2d) Process all small image blocks according to the processing procedure of image block A, and finally generate N. x *N y *N z N small image blocks, after being deformed, are placed in a folder. x N y N z This represents the number of blocks that can be taken along the X, Y, and Z axes.

[0070] Step (3) specifically includes the following steps:

[0071] (3a) First, all the small image patches of the target image form a coordinate system, ranging from [0,0,0 to [N]. x -1,N y -1,N z -1], therefore all small image blocks are divided into N z Layered data, each layer has N layers.x *N y There are several small image blocks, and let any layer k∈[0,N] be an example. z The Z-axis size of the small image patch [-1] is Z max The standard block size is set to [X]. b ,Y b Z b And satisfy X b =Y b =Z b =S per Under the condition of Z max The definition is as follows:

[0072]

[0073] Among them, Z n The Z-axis size of the target image;

[0074] (3b) Set a control variable Z during the splicing process. control Control variable Z control Used to control the amount of memory consumed when reading small image blocks;

[0075] (3c) When Z max ≤Z control At that time, read N of the kth layer x *N y Each small image block is read as Z. max For layer data, let the coordinates of any voxel within any small image block, i.e., image block A, be [x, y, z]. Then the voxel coordinates of this voxel in the target image are [x+X]. offset y+Y offset z+Z offset The offset of the origin of the starting coordinates of image block A in the target image is set to [X]. offset Y offset Z offset The image is stitched together using the X and Y axis offsets of each small image block, ultimately generating a Z-shaped image. max Zhang Image;

[0076] (3d) When Z max >Z control At this time, the processing steps of reading and stitching images K times are required, and their definitions are as follows:

[0077]

[0078] Each time, N is read from the k-th layer. x *N y Each small image block is read as Z. controlThe data from the layer is processed, and the image stitching process mentioned in step (3c) is performed to generate Z. control Zhang slices;

[0079] (3e) N z The layer data are processed in steps (3c) and (3d) to finally generate Z. n Zhang 2D continuous sequence images.

[0080] Step (2a) specifically includes the following steps:

[0081] (2a1) First, obtain the size of the target image. Since the ultra-large 3D image data is stored as a continuous 2D sequence of images in a folder and arranged in Z-axis order, read any image in the folder and obtain its X and Y axis dimensions, denoted as X. n Y n This refers to the X and Y axis dimensions of the target image, and the number of images read from the folder is denoted as Z. n This refers to the Z-axis size of the target image;

[0082] (2a2) Define a standard block size [X] b ,Y b Z b And satisfy X b =Y b =Z b =S per Let N x N y N z This represents the number of blocks that can be taken along the X, Y, and Z axes. These blocks include both standard and non-standard blocks, as defined below:

[0083]

[0084]

[0085]

[0086] (2a3) If the target image boundary is encountered, the size of the small image block is [X]. b′ ,Y b′ Z b′ The specific formula is as follows:

[0087] X b′ =X n -N x *X b

[0088] Y b′ =Y n -N y *Yb

[0089] Z b′ =Z n -N z *Z b

[0090] (2a4) The final number of small image blocks is determined to be N based on the target image size. sum Its definition is as follows:

[0091] N sum =N x *N y *N z .

[0092] In step (2b), the image deformation of each small image block specifically includes the following steps:

[0093] (2b1) After the target image is divided into blocks, a coordinate system is set according to the origin of the target image. Each small image block then obtains a corresponding image block number, i.e., the number range is:

[0094] [0,0,0] to [N] x -1,N y -1,N z -1]

[0095] For any image patch A, let its corresponding image patch number be [a,b,c], and its image patch size be [X]. A ,Y A Z A If ], then its X-axis dimension is X. b or X b′ Its Y-axis dimension is Y b Or Y b′ Its Z-axis dimension is Z b or Z b′ ;

[0096] (2b2) ​​Since each small image block is processed, the offset of the origin of the starting coordinates of image block A in the target image is set to [X]. offset Y offset Z offset ], which is defined as:

[0097] X offset =a*X b

[0098] Y offset =b*Y b

[0099] Z offset =c*Z b

[0100] The X and Y axis offsets of each image block are recorded and stored in two lists respectively;

[0101] (2b3) Based on the offset obtained in step (2b2), calculate the coordinates of the eight vertices of image block A, such as... Figure 2 As shown, for example, the coordinates of its vertices P1 and P7 can be represented as:

[0102]

[0103]

[0104] Among them, [X offset ,Y offset Z offset X is defined as the offset of the origin of the starting coordinates of image patch A in the target image. A Y A Z A Let X, Y, and Z be the dimensions of image patch A, and the same applies to other vertices; Let P1 be the coordinate of the vertex. Let P7 be the coordinates of vertex P7;

[0105] (2b4) Use the affine transformation matrix to convert the coordinates of the eight vertices of image block A into the coordinates of the eight vertices of the original image, and sort the coordinates in the three directions of X, Y and Z. Select the point with the largest and smallest values ​​in these three directions. The difference between the two values ​​is the range of the selected image block in the original image.

[0106] Step (2c) specifically includes the following steps:

[0107] (2c1) First, allocate N to the CUDA kernel function. Z ×N Y N threads, of which N Z N represents the Z-axis size of the 3D image to be registered. Y Set the Y-axis size of the 3D image to be registered, and then fix the blockDim in the kernel function to a 2D block with a size of 32*32;

[0108] (2c2) Each thread has a unique row number and column number, with the row number ranging from 0 to N. Z Its column number ranges from 0 to N. Y The size of the Z and Y axes corresponds to that of the three-dimensional image data;

[0109] (2c3) Then put the CUDA kernel function into a space of size 0 to N. X In the loop, N XThe X-axis size is the size of the 3D image to be registered. Finally, each thread processes the work of each voxel in the 3D image to be registered. Assuming that this is the kth loop, the thread with thread number (m, n) processes the voxel corresponding to the X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the 3D voxel. That is, each thread number corresponds to one voxel in the 3D image to be registered.

[0110] This invention utilizes the principle of segmenting the target image and leverages the parallel architecture of the GPU to accelerate the calculation process of linear image interpolation. This not only solves the problem of deforming ultra-large 3D images but also effectively improves the efficiency of image registration. Secondly, experiments have shown that this method can process ultra-large 3D image data larger than twice the computer memory on a platform with an i5-10400F CPU, an RTX-3060 GPU, and 16GB of memory, achieving an acceleration ratio of approximately 11 times. This significantly shortens the deformation registration time, as shown in Table 1, which illustrates the running time and processed image size results of this invention.

[0111] Table 1

[0112]

[0113] As can be seen from Table 1, based on a platform with an i5-10400F CPU, an RTX-3060 GPU, and 16GB of memory, the ultra-large 3D image size that can be processed is more than twice the size of the computer memory, and the speedup ratio can reach about 11 times, while only occupying 2.8GB of computer memory.

[0114] In summary, this invention utilizes the principle of segmenting the target image and leverages the parallel architecture of the GPU to accelerate the calculation process of linear image interpolation. This not only solves the problem of deforming ultra-large 3D images but also effectively improves the efficiency of image registration. Experiments show that on a platform with an i5-10400F CPU, an RTX-3060 GPU, and 16GB of memory, this method can process ultra-large 3D image data with more than twice the computer memory and achieves an acceleration ratio of approximately 11 times, greatly shortening the deformation registration time.

Claims

1. A method for deforming ultra-large 3D images based on CUDA acceleration, characterized in that: The method includes the following steps in sequence: (1) Calculate the affine transformation matrix of the two sets of registration feature points through affine transformation; (2) Divide the target image into blocks, and perform image deformation on each small image block obtained after the block division to obtain the deformed 3D image block: using the affine transformation matrix, the three-dimensional image voxel traversal parallelization operation is realized through CUDA, and trilinear interpolation is performed on each image block. (3) Perform image stitching operation for each small image block: read the deformed 3D image blocks and stitch them into a 2D continuous sequence of images of the original image size; Step (2) specifically includes the following steps: (2a) Divide the target image into blocks to obtain small image blocks; (2b) Perform deformation processing on each small image block to obtain the deformed 3D image block; (2c) Read the deformed 3D image block into memory, use CUDA to implement parallel operation of 3D image voxel traversal, use GPU-accelerated trilinear interpolation to process any small image block, i.e., image block A, and finally process the image block. The offset is used to name and save the file, and the file is saved in 3D image format; (2d) According to image blocks The processing steps involve processing all small image patches to ultimately generate... The image is divided into small image blocks after being deformed, and placed in a folder. This represents the number of blocks that can be taken along the X, Y, and Z axes; Step (3) specifically includes the following steps: (3a) First, all the small image patches of the target image form a coordinate system with a range of [missing information]. Therefore, all small image blocks are divided into Layered data, each layer has A small image block, with any layer The Z-axis size of the small image patch is The standard block size is set as follows: And satisfy Under the conditions, The definition is as follows: , in, The Z-axis size of the target image; (3b) Set a control variable during the splicing process. Control variables Used to control the amount of memory consumed when reading small image blocks; (3c) When At that time, read the first Layer Each small image block is read. The data of the layer, let the coordinates of any pixel within any small image block, i.e., image block A, be . Then the voxel coordinates of this voxel in the target image are Image blocks The offset of the origin of the starting coordinate system in the target image is set to The image is stitched together using the X and Y axis offsets of each small image block, ultimately generating... Zhang Image; (3d) When At that time, it is necessary to perform The image reading and stitching process is defined as follows: , Each time the first read Layer of common Each small image block is read. The data from the layer is processed, and the image stitching process mentioned in step (3c) is performed to generate... Zhang slices; (3e) will The layer data are processed in steps (3c) and (3d) to finally generate... Zhang 2D continuous sequence images.

2. The CUDA-accelerated ultra-large 3D image deformation method according to claim 1, characterized in that: Step (1) specifically refers to: obtaining two sets of data that each has... The registration feature points are used to normalize the two sets of registration feature points. Then, the affine transformation matrix is ​​calculated according to the affine transformation formula. Finally, the affine transformation matrix is ​​inversely normalized to obtain the final affine transformation matrix as follows: , in, ,..., All are affine parameters.

3. The CUDA-accelerated ultra-large 3D image deformation method according to claim 1, characterized in that: Step (2a) specifically includes the following steps: (2a1) First, obtain the size of the target image. Since the ultra-large 3D image data is stored as a continuous 2D sequence of images in a folder and ordered by the Z-axis, read any image from the folder and obtain its X and Y axis dimensions, denoted as... This refers to the X and Y axis dimensions of the target image, and the number of images read from the folder is recorded as follows. This refers to the Z-axis size of the target image; (2a2) Define a standard block size And satisfy ,set up The number of blocks that can be taken along the X, Y, and Z axes is defined as follows: , , , (2a3) If the target image boundary is encountered, then the size of the small image block is... The specific formula is as follows: , , , (2a4) The final number of small image blocks is determined based on the target image size. Its definition is as follows: 。 4. The CUDA-accelerated method for ultra-large 3D image deformation according to claim 1, characterized in that: In step (2b), the image deformation of each small image block specifically includes the following steps: (2b1) After the target image is divided into blocks, a coordinate system is set according to the origin of the target image. Each small image block then obtains a corresponding image block number, i.e., the number range is: , For any image patch Let its corresponding image block number be Let its image patch size be Then its X-axis dimension is or Its Y-axis size is or Its Z-axis dimension is or ; (2b2) ​​Since the processing is performed on each small image block, the image block The offset of the origin of the starting coordinate system in the target image is set to Its definition is: , , , The X and Y axis offsets of each image block are recorded and stored in two lists respectively; (2b3) Calculate the image patch based on the offset obtained in step (2b2). The coordinates of the eight vertices; (2b4) Using an affine transformation matrix to transform image patches The coordinates of the eight vertices are converted into the coordinates of the eight vertices of the original image. The coordinates are then sorted along the X, Y, and Z axes. The point with the largest and smallest coordinates along these three axes are selected, and the difference between the two is the range of the selected image patch in the original image.

5. The CUDA-accelerated method for ultra-large 3D image deformation according to claim 1, characterized in that: Step (2c) specifically includes the following steps: (2c1) First, allocate CUDA kernel functions. 1 thread, of which The Z-axis size of the 3D image to be registered. Set the Y-axis size of the 3D image to be registered, and then fix the blockDim in the kernel function to a 2D block with a size of 32*32; (2c2) Each thread has a unique row number and column number, with the row number ranging from 0 to 1. Its column number range is 0 to The size of the Z and Y axes corresponds to that of the three-dimensional image data; (2c3) Then put the CUDA kernel function into a space of size 0 to In the loop, The X-axis size is the size of the 3D image to be registered. Finally, each thread processes the work of each voxel in the 3D image to be registered. Assuming that this is the kth loop, the thread with thread number (m, n) processes the voxel with X-axis coordinate k, Y-axis coordinate n, and Z-axis coordinate m in the 3D voxel. That is, each thread number corresponds to one voxel in the 3D image to be registered.