CN101980304A

CN101980304A - Three-dimensional digital volume image distortion measuring method

Info

Publication number: CN101980304A
Application number: CN 201010520673
Authority: CN
Inventors: 黄建永; 潘晓畅; 李姗姗; 彭小玲; 熊春阳; 方竞
Original assignee: Peking University
Current assignee: Peking University
Priority date: 2010-10-20
Filing date: 2010-10-20
Publication date: 2011-02-23

Abstract

The invention relates to a three-dimensional digital volume image distortion measuring method, which comprises the following steps of: (1) selecting scattered sampling points on a digital volume image before distortion; (2) determining a reference sub-volume and the dimension of a searched area; (3) performing three-dimensional image-related operation on the reference sub-volume and the searched area based on a three-dimensional Fast Fourier Transform method; (4) establishing an image gray three-dimensional summing table and an image energy three-dimensional summing table; (5) solving a three-dimensional zero-mean normalization cross correlation coefficient matrix by using the operational result in the step (3) and the three-dimensional summing tables; (6) calculating subvoxel displacement values of the sampling points by using a gradient-based three-dimensional subvoxel displacement location algorithm; (7) calculating accurate displacement values at all the sampling points according to the steps (3) to (6) so as to obtain a three-dimensional displacement field of the whole distorted volume image relative to the volume image before distortion; and (8) calculating a three-dimensional strain field of the digital volume image according to the Lagrangian strain tensor theory in mechanics of continuous media. The three-dimensional distortion of the digital volume image can be accurately and efficiently measured.

Description

Three-dimensional digital volume image deformation measurement method

Technical Field

The invention relates to a measuring method, in particular to a deformation measuring method for quantitatively analyzing a three-dimensional digital volume image.

Background

The basic principle of the related art of two-dimensional digital images is to compare a digital image before deformation (reference image) with an image after deformation (deformation image), i.e. to compare the reference image with the deformation image. And then a series of cross correlation operations are carried out to obtain a displacement field and a strain field of the deformed image relative to the image before deformation. At present, the traditional two-dimensional digital image correlation method becomes a classical experimental mechanics technology which is widely applied to the fields of material mechanics behavior testing, micro-scale deformation measurement, micro-electro-mechanical system (MEMS) structure dynamic characterization, biological tissue mechanical property evaluation and the like. On the basis of the two-dimensional Digital image Correlation method, Bay and its co-workers proposed the concept of Digital Volume Correlation (DVC) in 1999 (Bay B.K., Smith T.S., Fyhrie D.P., and Saad M., Digital Volume Correlation: Three-dimensional mapping using X-ray tomogry, Experimental Mechanics, 1999, 39(3), 217-226; Bell, Smith, Share, Seawatt, Digital Volume image Correlation: Three-dimensional strain characterization using X-ray tomography, Experimental Mechanics, 1999, 39(3), 217-226), and was first applied to the deformation measurement of Three-dimensional X-ray microstructure in osteoporosis. The basic principle of the digital volume correlation method is similar to that of the digital image correlation method, and the core idea is to perform three-dimensional image cross-correlation operation between a specific Reference sub-volume (Reference Subvolume) in a deformed precursor image and all possible Target subvolumes (Target subvolumes) in a corresponding search area in the deformed volume image, and determine the most possible position of the Reference Subvolume in the deformed volume image through the extreme condition of a correlation function, so as to obtain a three-dimensional displacement field and a strain field of the digital volume image.

With the development of Micro-computed Tomography (μ CT), Magnetic Resonance Imaging (MRI), and other three-dimensional microstructure topography detection technologies in recent years, quantitative analysis of three-dimensional deformation of porous honeycomb materials (such as bones, wooden materials, clay materials, foam structure materials, and the like) by means of a digital volume correlation method is started. With the rapid development of laser-confocal microscopy (LSCM), research was undertaken to investigate three-dimensional deformation of soft materials with random markers and better light transmission (e.g., collagen materials with randomly distributed protein fibers) by using a digital volume-dependent method (Roeder b.a., Kokini k., Robinson j.p., and volume-Harbin s.l., Local, three-dimensional structure measures with transformed extracellular matrix structures, Journal of biological engineering, 2004, 126(6), 699-708; jord, cockayne, Robinson, waruck-708, strain measurement of extracellular matrix structures with Local three-dimensional large deformation, Journal of biological engineering, 2004, 69126 (6), 699). However, compared with the digital image correlation method which is developed more mature, the digital volume correlation method still needs to be developed and perfected in the aspects of solving precision and calculating efficiency, and the application field of the digital volume correlation method still needs to be further expanded.

At present, the digital volume correlation method mainly only considers the Voxel (also called Voxel) correlation matching operation between a reference sub-volume and a target sub-volume in the application process, and usually considers at most the rigid body translation and rotation effects between the sub-volumes in the matching process, and the sub-Voxel (subvolume) positioning algorithm framework of the digital volume correlation method is not established systematically, and particularly the research of the three-dimensional image sub-Voxel positioning algorithm based on laser confocal imaging is blank. On the other hand, as can be seen from the basic principle of the digital volume correlation method, for each selected reference sub-volume, three-dimensional volume correlation operations are required with all possible target sub-volumes within the corresponding search region one by one, so that it is not difficult to find that the correlation operations are much more computationally intensive than the search calculations of the digital image correlation method under similar circumstances. For example, when the sub-volume size is 32 × 32 × 32 voxels and the search region is 62 × 62 × 62 voxels, at least 30 × 30 × 30 ═ 27000 three-dimensional volume correlation operations are required between the reference sub-volume and the target sub-volume, and at least three-dimensional summation operations are required for each volume correlation calculation. Rough estimates indicate that the digital volume correlation computation in this case is approximately three orders of magnitude greater than the digital image correlation computation in a similar case (when the sub-region size is 32 x 32 pixels and the search region size is 62 x 62 pixels). The complexity of the above conventional digital volume correlation method in terms of calculation has become a major obstacle that restricts the further widespread application of the conventional digital volume correlation method, and particularly, the calculation efficiency of the conventional digital volume correlation method cannot meet the requirement in the case of large sample statistical analysis.

Disclosure of Invention

In view of the above problems, the present invention aims to provide a method for measuring the deformation of a three-dimensional digital volume image with high calculation efficiency and sub-voxel measurement accuracy.

In order to achieve the purpose, the invention adopts the following technical scheme: a three-dimensional digital volume image deformation measurement method comprises the following steps: (1) taking the three-dimensional digital volume image without deformation as a reference volume image F (x, y, z), and taking the three-dimensional digital volume image with deformation as a deformation volume image G (x, y, z); selecting a series of discrete sampling points on a reference volume image F (x, y, z), and setting the coordinate of any one sampling point as (alpha, beta, A); (2) respectively selecting L from the sampling points in the step (1)_x×L_y×L_zA rectangular parallelepiped of individual voxels as a reference subvolume f (x, y, z); in the deformed volume image G (x, y, z), M is selected from the sample point coordinates (alpha, beta, A) as the center_x×M_y×M_zA rectangular parallelepiped of each voxel is used as a search region g (x, y, z), and M_x＝M_y＝M_z＞L_x＝L_y＝L_z(ii) a (3) For each sample point (α, β, γ),and respectively carrying out three-dimensional image correlation operation on the corresponding reference sub-volume f (x, y, z) and the corresponding search area g (x, y, z), wherein the operation result P is as follows:

wherein represents a complex conjugate; FT₃[·]And

respectively representing three-dimensional fast Fourier forward transformation and inverse transformation;

and

representing a reference subvolume f (x, y, z) image grayscale three-dimensional matrix and a search area g (x, y, z) volume image grayscale three-dimensional matrix respectively; (4) respectively establishing a three-dimensional summation table S of the gray levels of the deformed digital volume images by utilizing a fast recursion method according to the image gray level value of the search area g (x, y, z) corresponding to each sampling point (alpha, beta, gamma)_gAnd image energy three-dimensional summation table(5) The three-dimensional summation table S established by the step (4)_gAndperforming fast table lookup operation, and calculating a three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w) ] between the reference subvolume f (x, y, z) corresponding to the sampling point (alpha, beta, gamma) and the corresponding search area g (x, y, z) according to the correlation operation result of the step (3)]_{(α，β，γ)}(ii) a (6) Performing sub-voxel interpolation operation near the position of the maximum peak of the zero-mean normalized cross-correlation coefficient matrix by adopting a gradient-based three-dimensional sub-voxel displacement positioning algorithm to obtain the accurate position coordinates (U, V, W) of the sampling points (alpha, beta, gamma) in the reference digital volume image in the corresponding deformed volume image, wherein U is U + delta U, V is V + delta V, W is W + delta W, and (u, v, w) represents a matrix of correlation coefficients [ C (u, v, w) normalized by the three-dimensional zero mean]_{(α，β，γ)}The integral voxel displacement value determined by the medium maximum element (delta u, delta v, delta w) is a sub-voxel displacement value calculated by a three-dimensional sub-voxel displacement positioning algorithm based on gradient; (7) repeating the steps (3) to (6), and calculating the accurate positions of the sampling points in all the reference volume images in the deformed digital volume image so as to obtain a three-dimensional displacement field of the whole deformed volume image relative to the reference volume image; (8) calculating a three-dimensional strain field of the digital volume image according to a Lagrange strain tensor theory in continuous medium mechanics:wherein, i is more than or equal to 1, j, m is less than or equal to 3, U₁＝U，U₂＝V，U₃＝W。

In the step (2), in the step L_x×L_y×L_zIn a rectangular parallelepiped of individual voxels, L_x＝L_y＝L_zAnd the integers are integers of 30-60, and the unit is a voxel; at the M_x×M_y×M_zIn a rectangular parallelepiped of individual voxels, M_x＝M_y＝M_zAnd they take integers between 35 and 100, the unit is voxel, and M_x＝M_y＝M_z＞L_x＝L_y＝L_z。

In the step (4), the three-dimensional summation table S_gAnd

respectively:

s_g(x，y，z)＝g(x，y，z)+s_g(x，y-1，z)+s_g(x-1，y，z)+s_g(x-1，y-1，z-1)

+s_g(x，y，z-1)-s_g(x，y-1，z-1)-s_g(x-1，y，z-1)-s_g(x-1，y-1，z)’

s_{g^{2}} (x, y, z) = {[g (x, y, x)]}^{2} + s_{g^{2}} (x, y - 1, z) + s_{g^{2}} (x - 1, y, z) + s_{g^{2}} (x - 1, y - 1, z - 1)

+ s_{g^{2}} (x, y, z - 1) - s_{g^{2}} (x, y - 1, z - 1) - s_{g^{2}} (x - 1, y, z - 1) - s_{g^{2}} (x - 1, y - 1, z),

wherein, the three-dimensional summation table S_gAnd

respectively comprise M_x×M_y×M_zAn element; g (x, y, z) represents the gray value corresponding to the voxel with the central coordinate (x, y, z) in the search area; and x, y, z are integers, and when x, y, z is less than or equal to 0,

S_{g} = S_{g^{2}} = 0 .

in the step (5), the three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w)]_{(α，β，γ)}Comprises the following steps:

wherein,

<math><mrow><mi>Q</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo>)</mo></mrow><mo>=</mo><mover><mi>f</mi><mo>&OverBar;</mo></mover><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>,</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><msub><mi>L</mi><mi>y</mi></msub><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><mo>[</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>]</mo><mo>×</mo><mo>[</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>u</mi></mrow><mrow><msub><mi>L</mi><mi>x</mi></msub><mo>+</mo><mi>u</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>v</mi></mrow><mrow><msub><mi>L</mi><mi>y</mi></msub><mo>+</mo><mi>v</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>w</mi></mrow><mrow><msub><mi>L</mi><mi>z</mi></msub><mo>+</mo><mi>w</mi></mrow></munderover><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>]</mo><mo>,</mo></mrow></math>

<math><mrow><mi>G</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo>)</mo></mrow><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><msup><mrow><mo>[</mo><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>,</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow><mo>-</mo><msub><mover><mi>g</mi><mo>&OverBar;</mo></mover><mi>uvw</mi></msub><mo>]</mo></mrow><mn>2</mn></msup><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>u</mi></mrow><mrow><msub><mi>L</mi><mi>x</mi></msub><mo>+</mo><mi>u</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>v</mi></mrow><mrow><msub><mi>L</mi><mi>y</mi></msub><mo>+</mo><mi>v</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>w</mi></mrow><mrow><msub><mi>L</mi><mi>z</mi></msub><mo>+</mo><mi>w</mi></mrow></munderover><msup><mi>g</mi><mn>2</mn></msup><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><msub><mi>L</mi><mi>y</mi></msub><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><msup><mrow><mo>[</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>u</mi></mrow><mrow><msub><mi>L</mi><mi>x</mi></msub><mo>+</mo><mi>u</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>v</mi></mrow><mrow><msub><mi>L</mi><mi>y</mi></msub><mo>+</mo><mi>v</mi></mrow></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn><mo>+</mo><mi>w</mi></mrow><mrow><msub><mi>L</mi><mi>z</mi></msub><mo>+</mo><mi>w</mi></mrow></munderover><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>,</mo></mrow></math>

<math><mrow><mi>F</mi><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mrow><mo>[</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mover><mi>f</mi><mo>&OverBar;</mo></mover><msup><mo>]</mo><mn>2</mn></msup><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><msup><mi>f</mi><mn>2</mn></msup><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><msub><mi>L</mi><mi>y</mi></msub><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><msup><mrow><mo>[</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>,</mo></mrow></mrow></math>

and is

<math><mrow><mover><mi>f</mi><mo>&OverBar;</mo></mover><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><mo>×</mo><msub><mi>L</mi><mi>y</mi></msub><mo>×</mo><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>;</mo></mrow></math>

<math><mrow><msub><mover><mi>g</mi><mo>&OverBar;</mo></mover><mi>uvw</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><mo>×</mo><msub><mi>L</mi><mi>y</mi></msub><mo>×</mo><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>,</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow><mo>;</mo></mrow></math>

P (u, v, w) is rapidly solved by the three-dimensional fast Fourier transform method in the step (3), i.e.Solving the triple summation items in Q (u, v, w) and G (u, v, w) by a quick table look-up mode in the step (5); f is solved by arithmetic operations.

Due to the adoption of the technical scheme, the invention has the following advantages: 1. according to the invention, a three-dimensional fast Fourier transform method is adopted to carry out fast cross correlation operation on the reference sub-volumes and the corresponding search areas, a three-dimensional image gray level summation table and a three-dimensional image energy summation table are established in a fast recursion mode, and then a table look-up mode is used to quickly calculate the triple summation expression, so that efficient solution of the three-dimensional zero-mean normalized cross correlation coefficient is realized, and further the integer element displacement result of the sampling point is quickly obtained. 2. The invention adopts a fast table look-up mode to calculate the triple summation expression in the three-dimensional zero mean value normalization cross correlation coefficient, and introduces the three-dimensional fast Fourier transform technology to carry out correlation operation, thereby effectively reducing the complexity of the displacement calculation of the sampling point, and particularly improving the solving efficiency of the digital volume image correlation operation under the condition of carrying out large-scale volume correlation operation or having harsh requirements on the calculation time. 3. According to the invention, the three-dimensional sub-voxel positioning algorithm based on the spatial gradient is adopted to carry out sub-voxel displacement solving on the digital volume image, so that the solving precision of the displacement field can reach the sub-voxel magnitude, and the sub-voxel positioning algorithm does not need to carry out complicated iterative calculation, so that the calculation efficiency is higher. The invention can be widely applied to the deformation measurement analysis of the three-dimensional digital volume image.

Drawings

FIG. 1 is a schematic overall flow diagram of the present invention;

FIG. 2 is a schematic diagram of the digital volume image correlation principle of the present invention;

FIG. 3 is a schematic diagram of three-dimensional digital volume image effects before and after deformation after undergoing vertical single-cycle compression according to an embodiment of the present invention;

FIGS. 4 a-4 c are schematic diagrams of full-field three-dimensional displacement field effects calculated according to a first embodiment of the present invention;

fig. 5a to 5b are schematic diagrams illustrating comparison of three-dimensional deformation effect of agarose-fluorescence digital volume images captured by a confocal laser microscope and measured by the present invention according to a first embodiment of the present invention, wherein the units are voxels, fig. 5a is a three-dimensional agarose-fluorescence substrate digital volume deformation image captured by a confocal laser microscope, and fig. 5b is a displacement field calculated by the three-dimensional digital volume image deformation measuring method according to the present invention.

Detailed Description

The invention is described in detail below with reference to the figures and examples.

As shown in fig. 1 and fig. 2, the three-dimensional digital volume image deformation measurement method of the present invention is based on the inherent characteristics of digital volume image three-dimensional deformation search calculation, and by introducing a three-dimensional fast fourier transform method and combining a three-dimensional summation table method, the steps are as follows:

1) taking the three-dimensional digital volume image without deformation as a reference volume image F (x, y, z), and taking the three-dimensional digital volume image with deformation as a deformation volume image G (x, y, z); selecting a series of discrete sampling points in a reference volume image F (x, y, z) according to a certain rule (for example, selecting grid intersection points as sampling points by a method of dividing a stereo grid), and setting any one of the sampling points to be represented by coordinates (alpha, beta, gamma);

2) respectively selecting a sample containing L by taking a series of sampling points in the step 1) as a center_x×L_y×L_zA rectangular parallelepiped of individual voxels serves as a reference subvolume f (x, y, z), wherein L is usually given_x＝L_y＝L_zAnd they generally take integers between 30 and 60, and the unit is voxel; similarly, in the deformed digital volume image (also called deformed volume image) G (x, y, z), the coordinates (α, β, γ) of the sampling points previously selected in the reference volume image are respectively selected to include M_x×M_y×M_zA rectangular parallelepiped of individual voxels serves as the search region g (x, y, z), where M is usually assumed_x＝M_y＝M_zAnd they generally take integers between 35 and 100 in units of voxels, and M_x＝M_y＝M_z＞L_x＝L_y＝L_z(as shown in FIG. 2);

3) for each sampling point (α, β, γ), performing three-dimensional image correlation operation on the corresponding reference sub-volume f (x, y, z) and the corresponding search area g (x, y, z), wherein an operation result P is:

P = {FT}_{3}^{- 1} [{FT}_{3} (f) {FT}_{3}^{*} (g)], - - - (1)

wherein represents a complex conjugate; FT₃[·]And

and

representing a reference subvolume f (x, y, z) image grayscale three-dimensional matrix and a search area g (x, y, z) volume image grayscale three-dimensional matrix respectively;

the present invention uses the above-described three-dimensional fast fourier transform method to perform three-dimensional image correlation operations on the reference subvolume f (x, y, z) and the search area g (x, y, z), which can reduce the computational complexity of the prior art from 0[ (L)_x×L_y×L_z)×(M_x-L_x+1)×(M_y-L_y+1)×(M_z-L_z+1)]Reduced to O [ (M)_x×M_y×M_z)×log₂(M_x×M_y×M_z)]；

4) Respectively establishing a three-dimensional summation table S of the gray levels of the deformed digital volume images by utilizing a fast recursion method according to the image gray level value of the search area g (x, y, z) corresponding to each sampling point (alpha, beta, gamma)_gAnd image energy three-dimensional summation table

Namely:

(2)

+s_g(x，y，z-1)-s_g(x，y-1，z-1)-s_g(x-1，y，z-1)-s_g(x-1，y-1，z)，

s_{g^{2}} (x, y, z) = {[g (x, y, z)]}^{2} + s_{g^{2}} (x, y - 1, z) + s_{g^{2}} (x - 1, y, z) + s_{g^{2}} (x - 1, y - 1, z - 1),

(3)

+ s_{g^{2}} (x, y, z - 1) - s_{g^{2}} (x, y - 1, z - 1) - s_{g^{2}} (x - 1, y, z - 1) - s_{g^{2}} (x - 1, y - 1, z)

in the formula, a three-dimensional summation table S_gAndrespectively comprise M_x×M_y×M_zAn element; g (x, y, z) here represents the gray value corresponding to the voxel with central coordinate (x, y, z) in the search area; x, y and z are integers, and when x, y and z are less than or equal to 0,

S_{g} = S_{g^{2}} = 0;

5) the three-dimensional summation table S established by the step 4)_gAnd

by performing a fast table lookup operation, the following triple summation expression can be obtained quickly, that is:

= s_{g} (L_{x} + u, L_{y} + v, L_{z} + w) + s_{g} (L_{x} + u, v, w) + s_{g} (u, L_{y} + v, w) + s_{g} (u, v, L_{z} + w), - - - (4)

- s_{g} (L_{x} + u, v, L_{z} + w) - s_{g} ({u, L}_{y} + v, L_{z} + w) - s_{g} (L_{x} + u, L_{y} + v, w) - s_{g} (u, v, w)

= s_{g^{2}} (L_{x} + u, L_{y} + v, L_{z} + w) + s_{g^{2}} (L_{x} + u, v, w) + s_{g^{2}} (u, L_{y} + v, w) + s_{g^{2}} (u, v, L_{z} + w), - - - (5)

- s_{g^{2}} (L_{x} + u, v, L_{z} + w) - s_{g^{2}} ({u, L}_{y} + v, L_{z} + w) - s_{g^{2}} (L_{x} + u, L_{y} + v, w) - s_{g^{2}} (u, v, w)

accordingly, according to the correlation operation result in the step 3), calculating a three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w) ] between the reference sub-volume f (x, y, z) corresponding to the sampling point (alpha, beta, gamma) and the corresponding search area g (x, y, z)]_{(α，β，γ)}(ii) a The position of the maximum peak value of the matrix is the integral voxel displacement of the current sampling point (alpha, beta, gamma);

6) performing gradient-based three-dimensional sub-voxel displacement positioning algorithm near the position of the maximum peak value of the zero-mean normalized cross-correlation coefficient matrixPerforming sub-voxel interpolation operation to obtain accurate position coordinates (U, V, W) of sampling points (alpha, beta, gamma) in the reference digital volume image in the volume image after corresponding deformation; where U + Δ U, V + Δ V, W + W, and (U, V, W) represent a matrix of correlation coefficients normalized by a three-dimensional zero mean value [ C (U, V, W)]_(αβγ)The global voxel displacement value determined by the medium maximum element, (Δ u, Δ v, Δ w) is a sub-voxel displacement value calculated by a gradient-based three-dimensional sub-voxel displacement localization algorithm, and can be directly calculated by the following analytical formula:

in the formula, g_x、g_vAnd g_zCan be solved by the gray-scale weighted calculation of the adjacent voxels of the digital volume image respectively, namely:

\{\begin{matrix} g_{x} = [g (x - 2, y, z) - 8 g (x - 1, y, z) + 8 g (x + 1, y, z) - g (x + 2, y, z)] / 12 \\ g_{y} = [g (x, y - 2, z) - 8 g (x, y - 1, z) + 8 g (x, y + 1, z) - g (x, y + 2, z)] / 12 \\ g_{z} = [g (x, y, z - 2) - 8 g (x, y, z - 1) + 8 g (x, y, z + 1) - g (x, y, z + 2)] / 12 \end{matrix}; - - - (7)

7) repeating the steps 3) to 6), and calculating the accurate positions of the sampling points in all the reference volume images in the deformed digital volume image so as to obtain a three-dimensional displacement field of the whole deformed volume image relative to the reference volume image;

8) calculating a three-dimensional strain field of the digital volume image according to a Lagrange strain tensor theory in continuous medium mechanics:

<math><mrow><msub><mi>ϵ</mi><mi>ij</mi></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>[</mo><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>i</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>j</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>+</mo><munderover><mi>Σ</mi><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mn>3</mn></munderover><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>m</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>m</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>]</mo><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow></math>

wherein, i is more than or equal to 1, j, m is less than or equal to 3, U₁＝U，U₂＝V，U₃＝W。

In the above step 5), the three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w)]_{(α，β，γ)}Obtained by the following steps:

the three-dimensional zero-mean normalized cross-correlation coefficient matrix between the reference subvolume f (x, y, z) corresponding to the sampling point (α, β, γ) and the corresponding search area g (x, y, z) has the basic expression:

<math><mrow><msub><mrow><mo>[</mo><mi>C</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>,</mo><mi>w</mi><mo>)</mo></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mi>α</mi><mo>,</mo><mi>β</mi><mo>,</mo><mi>γ</mi><mo>)</mo></mrow></msub><mo>=</mo><mfrac><mrow><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mo>[</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mover><mi>f</mi><mo>&OverBar;</mo></mover><mo>]</mo><mo>[</mo><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>,</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow><mo>-</mo><msub><mover><mi>g</mi><mo>&OverBar;</mo></mover><mi>uvw</mi></msub><mo>]</mo></mrow><mrow><msqrt><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><msup><mrow><mo>[</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mover><mi>f</mi><mo>&OverBar;</mo></mover><mo>]</mo></mrow><mn>2</mn></msup></msqrt><msqrt><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><msup><mrow><mo>[</mo><mi>g</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>,</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow><mo>-</mo><msub><mover><mi>g</mi><mo>&OverBar;</mo></mover><mi>uvw</mi></msub><mo>]</mo></mrow><mn>2</mn></msup></msqrt></mrow></mfrac><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow></math>

wherein:

order to

<math><mrow><mi>F</mi><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mrow><msup><mrow><mo>[</mo><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mover><mi>f</mi><mo>&OverBar;</mo></mover><mo>]</mo></mrow><mn>2</mn></msup><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><msup><mi>f</mi><mn>2</mn></msup><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>-</mo><mfrac><mn>1</mn><mrow><msub><mi>L</mi><mi>x</mi></msub><msub><mi>L</mi><mi>y</mi></msub><msub><mi>L</mi><mi>z</mi></msub></mrow></mfrac><msup><mrow><mo>[</mo><munderover><mi>Σ</mi><mrow><mi>x</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>x</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>y</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>y</mi></msub></munderover><munderover><mi>Σ</mi><mrow><mi>z</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>L</mi><mi>z</mi></msub></munderover><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>,</mo></mrow></mrow></math>

Equation (9) can be simplified to:

secondly, the expression P (u, v, w) in the step (i) can be directly solved by the three-dimensional fast Fourier transform method in the step (3), namelyThe triple summation items in the expressions Q (u, v, w) and G (u, v, w) can be solved by the fast table look-up method in the step 5); in addition, the above expression F can be solved by direct arithmetic operation, and since the expression F only needs to be calculated once for any sampling point (α, β, γ), the calculation time cost can be approximately ignored.

The measuring method of the present invention is further described below by way of specific examples.

The first embodiment is as follows: as shown in fig. 3, it is assumed that the vertical single-cycle uniform compressive strain previously applied by human is 2.0%, 5.0%, and 7.0%, respectively. As shown in fig. 4a to 4c, fig. 5a and 5b, the full-field three-dimensional displacement field calculated by the method for measuring the deformation of the three-dimensional digital volume image of the present invention, the three-dimensional agarose-fluorescence-based digital volume deformation image photographed by the confocal laser microscope, and the displacement field calculated by the method for measuring the deformation of the three-dimensional digital volume image of the present invention are shown. The error between the artificially pre-strained and the actual calculated strain is shown in table 1.

TABLE 1 artificial prestrain field and three-dimensional digital volume image deformation measuring method of the invention

Calculating to obtain the average error and standard deviation between strain fields

After the three-dimensional summation table method is adopted, the acceleration ratio (i.e. the ratio of the calculation efficiency) of the triple summation operation in the three-dimensional zero-mean normalized cross-correlation coefficient is obtained as shown in table 2.

TABLE 2 acceleration ratio of triple summation operation after the use of a summation table

Therefore, the three-dimensional digital volume image deformation measurement method can efficiently carry out digital volume image correlation operation, and accurately acquire the three-dimensional deformation field of the digital volume image, and mainly comprises a three-dimensional displacement field and a three-dimensional strain field.

The above embodiments are only for illustrating the present invention, and the scope of the present invention is not limited thereto, and any modifications and equivalent changes of the individual steps according to the principle of the present invention should not be excluded from the scope of the present invention based on the technical method of the present invention.

Claims

1. A three-dimensional digital volume image deformation measurement method comprises the following steps:

(1) taking the three-dimensional digital volume image without deformation as a reference volume image F (x, y, z), and taking the three-dimensional digital volume image with deformation as a deformation volume image G (x, y, z); selecting a series of discrete sampling points on a reference volume image F (x, y, z), and setting the coordinates of any one sampling point as (alpha, beta, gamma);

(2) respectively selecting L from the sampling points in the step (1)_x×L_y×L_zA rectangular parallelepiped of individual voxels as a reference subvolume f (x, y, z); in the deformed volume image G (x, y, z), M is selected from the sample point coordinates (alpha, beta, gamma) as the center_x×M_y×M_zA rectangular parallelepiped of each voxel is used as a search region g (x, y, z), and M_x＝M_y＝M_z＞L_x＝L_y＝L_z；

(3) For each sampling point (α, β, γ), performing three-dimensional image correlation operation on the corresponding reference sub-volume f (x, y, z) and the corresponding search area g (x, y, z), wherein an operation result P is:

P = {FT}_{3}^{- 1} [{FT}_{3} (f) {FT}_{3}^{*} (g)],

wherein represents a complex conjugate; FT₃[·]Andrespectively representing three-dimensional fast Fourier forward transformation and inverse transformation;

andrepresenting a reference subvolume f (x, y, z) image grayscale three-dimensional matrix and a search area g (x, y, z) volume image grayscale three-dimensional matrix respectively;

(4) respectively establishing a three-dimensional summation table S of the gray levels of the deformed digital volume images by utilizing a fast recursion method according to the image gray level value of the search area g (x, y, z) corresponding to each sampling point (alpha, beta, gamma)_gAnd image energy three-dimensional summation table

(5) The three-dimensional summation table S established by the step (4)_gAnd

performing fast table lookup operation, and calculating a three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w) ] between the reference subvolume f (x, y, z) corresponding to the sampling point (alpha, beta, gamma) and the corresponding search area g (x, y, z) according to the correlation operation result of the step (3)]_{(α，β，γ)}；

(6) Performing sub-voxel interpolation operation near the position of the maximum peak of the zero-mean normalized cross-correlation coefficient matrix by adopting a gradient-based three-dimensional sub-voxel displacement positioning algorithm to obtain the accurate position coordinates (U, V, W) of the sampling points (alpha, beta, gamma) in the reference digital volume image in the corresponding deformed volume image, wherein U is U + delta U, V is V + delta V, W is W + delta W, and U, V, W represent the three-dimensional zero-mean normalized correlation coefficient matrix [ C (U, V, W)]_{(α，β，γ)}The integral voxel displacement value determined by the medium maximum element (delta u, delta v, delta w) is a sub-voxel displacement value calculated by a three-dimensional sub-voxel displacement positioning algorithm based on gradient;

(7) repeating the steps (3) to (6), and calculating the accurate positions of the sampling points in all the reference volume images in the deformed digital volume image so as to obtain a three-dimensional displacement field of the whole deformed volume image relative to the reference volume image;

(8) calculating a three-dimensional strain field of the digital volume image according to a Lagrange strain tensor theory in continuous medium mechanics:

<math><mrow><msub><mi>ϵ</mi><mi>ij</mi></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>[</mo><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>i</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>j</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>+</mo><munderover><mi>Σ</mi><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mn>3</mn></munderover><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>m</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mfrac><mrow><mo>&PartialD;</mo><msub><mi>U</mi><mi>m</mi></msub></mrow><mrow><mo>&PartialD;</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>]</mo><mo>,</mo></mrow></math>

2. A method of deformation measurement of a three-dimensional digital volumetric image as defined in claim 1, characterized by: in the step (2), in the step L_x×L_y×L_zIn a rectangular parallelepiped of individual voxels, L_x＝L_y＝L_zAnd the integers are integers of 30-60, and the unit is a voxel; at the M_x×M_Y×M_zIn a rectangular parallelepiped of individual voxels, M_x＝M_v＝M_zAnd they take integers between 35 and 100, the unit is voxel, and M_x＝M_y＝M_z＞L_x＝L_y＝L_z。

3. A method of deformation measurement of a three-dimensional digital volumetric image as defined in claim 1, characterized by: in the step (4), the three-dimensional summation table S_gAnd

respectively:

s_g(x，y，z)＝g(x，y，z)+s_g(x，y-1，z)+s_g(x-1，y，z)+s_g(x-1，y-1，z-1)+s_g(x，y，z-1)-s_g(x，y-1，z-1)-s_g(x-1，y，z-1)-s_g(x-1，y-1，z)，

s_{g^{2}} (x, y, z) = {[g (x, y, x)]}^{2} + s_{g^{2}} (x, y - 1, z) + s_{g^{2}} (x - 1, y, z) + s_{g^{2}} (x - 1, y - 1, z - 1)

+ s_{g^{2}} (x, y, z - 1) - s_{g^{2}} (x, y - 1, z - 1) - s_{g^{2}} (x - 1, y, z - 1) - s_{g^{2}} (x - 1, y - 1, z),

wherein, the three-dimensional summation table S_gAnd

S_{g} = S_{g^{2}} = 0 .

4. a method of deformation measurement of a three-dimensional digital volumetric image as defined in claim 2, characterized in that: in the step (4), the three-dimensional summation tables Sg and Sg

Respectively:

s_g(x，y，z)＝g(x，y，z)+s_g(x，y-1，z)+s_g(x-1，y，z)+s_g(x-1，y-1，z-1)+s_g(x，y，z-1)-s_g(x，y-1，z-1)-s_g(x-1，y，z-1)-s_g(x-1，y-1，z)’

s_{g^{2}} (x, y, z) = {[g (x, y, x)]}^{2} + s_{g^{2}} (x, y - 1, z) + s_{g^{2}} (x - 1, y, z) + s_{g^{2}} (x - 1, y - 1, z - 1)

+ s_{g^{2}} (x, y, z - 1) - s_{g^{2}} (x, y - 1, z - 1) - s_{g^{2}} (x - 1, y, z - 1) - s_{g^{2}} (x - 1, y - 1, z),

wherein, the three-dimensional summation table S_gAnd

S_{g} = S_{g^{2}} = 0 .

5. a method of deformation measurement of a three-dimensional digital volume image as claimed in claim 1 or 2 or 3 or 4, characterized by: in the step (5), the three-dimensional zero-mean normalized cross-correlation coefficient matrix [ C (u, v, w)_{(α，β，γ)}Comprises the following steps:

wherein,

and is

P (u, v, w) is rapidly solved by the three-dimensional fast Fourier transform method in the step (3), i.e.

Solving the triple summation items in Q (u, v, w) and G (u, v, w) by a quick table look-up mode in the step (5); f is solved by arithmetic operations.