Next Article in Journal
A Review of Data Analytic Applications in Road Traffic Safety. Part 2: Prescriptive Modeling
Next Article in Special Issue
Improved Depth-of-Field Photoacoustic Microscopy with a Multifocal Point Transducer for Biomedical Imaging
Previous Article in Journal
SNPL: One Scheme of Securing Nodes in IoT Perception Layer
Previous Article in Special Issue
Active 3D Imaging of Vegetation Based on Multi-Wavelength Fluorescence LiDAR
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

3D Measurement of Human Chest and Abdomen Surface Based on 3D Fourier Transform and Time Phase Unwrapping

The Higher Educational Key Laboratory for Measuring & Control Technology and Instrumentations of Heilongjiang Province, Harbin University of Science and Technology, Harbin 150080, China
*
Author to whom correspondence should be addressed.
Sensors 2020, 20(4), 1091; https://doi.org/10.3390/s20041091
Submission received: 13 January 2020 / Revised: 11 February 2020 / Accepted: 16 February 2020 / Published: 17 February 2020
(This article belongs to the Special Issue Imaging Sensors and Applications)

Abstract

:
Monitoring respiratory movements is an effective way to improve radiotherapy treatments of thoracic and abdominal tumors, but the current approach is limited to measuring specific points in the chest and abdomen. In this paper, a dynamic three-dimensional (3D) measurement approach of the human chest and abdomen surface is proposed, which can infer tumor movement more accurately, so the radiotherapy damage to the human body can be reduced. Firstly, color stripe patterns in the RGB color model are projected, then after color correction, the collected stripe image sequences are separated into the three RGB primary color stripe image sequences. Secondly, a fringe projection approach is used to extract the folded phase combined 3D Fourier transform with 3D Gaussian filtering. By the relationship between adjacent fringe images in the time sequence, Gaussian filter parameters with individual characteristics are designed and optimized to improve the accuracy of wrapped phase extraction. In addition, based on the difference between the fractional parts of the folded phase error, one remainder equation can be determined, which is used for time-phase unwrapping. The simulation model and human experiments show that the proposed approach can obtain the 3D image sequences of the chest and abdomen surface in respiratory motion effectively and accurately with strong anti-interference ability.

1. Introduction

During radiotherapy, respiratory movements can cause tumors and normal tissues of the chest and abdomen to move at a certain frequency and amplitude. Sometimes, respiratory movements may affect the radiotherapy effect and even cause radiotherapy damage to the human body. In order to solve this respiratory motion problem, now the most effective real-time tracking method is to monitor extracorporeal respiratory movement. Based on extracorporeal respiratory movement, the respiratory movement of the tumor can be inferred, and then, the relative position of the target area and the field can be controlled by the radiotherapy system [1].
The fact that the respiratory movement of the tumor can be deduced from extracorporeal respiratory movement has been proven to be effective [2]. Based on this premise, dynamic 3D measurement of the chest and abdomen surface can be used to infer tumor movement more accurately [3].
Currently, optical measurement has become the most practical method to solve the problem of dynamic 3D measurement of the chest and abdomen surface. The optical method contains the following three types, point imaging, line imaging, and surface imaging, of which the surface imaging method is the best choice for dynamic 3D measurement. One image or many images can be captured by the surface imaging method. Single image acquisition methods include binocular vision [4], spatially encoded light [5], Fourier transform profile [6], etc. Multiple image acquisition methods include the phase-shift profile [7], modulation measurement profile [8], etc. Due to the fact that the multiple image acquisition method needs to capture many images, it has low efficiency and is not suitable for dynamic measurements. Thus, the one image acquisition method is a good solution for dynamic three-dimensional measurement.
Among these single image acquisition methods, the binocular vision method needs complex stereo matching and has low accuracy [9]; and the spatially encoded light method needs to be coded and decoded by the neighboring pixels, which makes the measurement resolution limited and may cause measurement failure in the case of surface height jump or shade [10]. Fluoroscopic real-time tumor tracking radiotherapy following 4D treatment planning was developed and shown to be feasible to improve the accuracy of the radiotherapy for mobile tumors [11]. In radiation therapy, the projection patterns need to be simple and continuous changes. Because the Fourier transform profile method has significant advantages in noise suppression and full field measurement and the sinusoidal fringe patterns projected by this method meet the requirements of radiotherapy, the Fourier transform profile method is suitable for measuring the 3D motion of the thoracic and abdominal surfaces. However, when only one stripe pattern is projected to measure the entire chest and abdomen surface, the difference in light intensity between adjacent pixels is small due to the large measurement range, and the anti-interference ability is low [12].
In this paper, a dynamic 3D measurement approach of the human chest and abdomen surface during respiration is proposed, which provides a basis for inferring and tracking tumor respiration movement during radiotherapy. This approach adopts a single color stripe pattern with three periods. Through combining one coded pattern with the three RGB primary colors, the sinusoidal stripe pattern with three different periods can be formed. During measurement, the projection pattern does not change, and the deformed stripe image of the chest and abdomen surface is collected in real time. Then, after color coupling correction and color separation, the single color deformation fringe images with three different periods can be formed. The proposed approach can obtain three deformed fringe images by one unchanged pattern. Taking each image sequence of single color fringe as a whole, the 3D Fourier fringe analysis (3D-FFA) method is used to extract the folded phase. This method has higher anti-interference ability. The three-frequency time phase unwrapping method is adopted. The absolute phase is obtained by the folding phase of three monochromatic fringes. This method has a large unwrapping range and strong anti-interference ability. According to the principle of triangulation, the 3D coordinates of the chest and abdomen surface are obtained from the absolute phase [13].
Section 2 introduces the proposed 3D measurement system. Section 3 describes the folding phase extraction method. Section 4 derives a three-frequency time phase unwrapping method. Section 5 shows and analyzes the experimental results, and the conclusions are given in Section 6.

2. 3D Measurement System Description

Figure 1 is a schematic diagram of the 3D measurement system for the human chest and abdomen surface. The method mainly includes the following five parts.
(1)
Pattern projection and image acquisition. The computer generates different periods of three RGB primary color cosine stripe patterns, and these three parts are combined into a composite color stripe pattern, then this pattern is projected to the chest and abdomen surface of the human. The camera captures stripe images of the chest and abdomen surface, which change with breathing movements at regular intervals to get a composite color stripe image sequence. The projection pattern in the proposed approach does not change, which can reduce the time of projection pattern conversion and setup. In addition, the measurement system only collects one composite color stripe pattern, which can decrease the image acquisition time. All these advantages can lay the foundation of the dynamic 3D measurement for the human chest and abdomen surface.
(2)
Image color correction and separation. For compound color stripe image sequences, color coupling correction and color separation should be made based on the correction matrix of each pixel, and the three RGB primary color fringe image sequences of different periods can be separated. Because the color coupling phenomenon exists at the coincident intersection of the three color channel spectral response curves in the 3CCD industrial camera [14], the color calibration based on hardware equipment should be completed before the measurement. That is to say, the projector projects four patterns of full red, full green, full blue, and full black to the chest and abdomen surface, and the four images will be captured by the camera. Using these four images, the correction matrix of each pixel is obtained according to the Casti illumination model [15].
(3)
Folded phase extraction. For each image in the RGB stripe image sequence, three-dimensional Fourier fringe analysis (3D-FFA) is used to extract the folding phase of each pixel and get the folded phase map of each image, and then, RGB folding phase map sequences can be formed. Fourier fringe analysis (FFA) is extended from one-dimensional Fourier fringe analysis (1D-FFA) to two-dimensional Fourier fringe analysis (2D-FFA) by using the properties of 2D fringe images. After this process, useful signals and interference can be separated better. This method becomes an effective measurement for 3D measurement of flat surfaces [16]. In this paper, the image sequences of the chest and abdomen surface are taken as a 3D one, which is analyzed by 3D Fourier transformation. Useful signals and interference can be separated further by increasing the time dimension, so as to reduce the influence of interference and improve the accuracy of measurement.
(4)
Folded phase unwrapping. According to the RGB folding phase diagram at the same time, with the proposed method of three-frequency time phase unwrapping in this paper, the folded phase is expanded into a continuous absolute phase, and the absolute phase diagram at that moment is obtained; thus, an absolute phase sequence can be formed. In the phase unwrapping method of this paper, the unwrapping operation depends on the difference of the decimal part of the measured folded phase, which can ensure that the absolute phase error does not exceed the folded phase error under certain conditions. In addition, we can judge whether there is any big error based on the absolute phase value, which can eliminate or reduce the effect of large absolute phase error by eliminating or interpolating operation. The phase unwrapping is achieved by solving the remainder equation set in the maximum range.
(5)
Three-dimensional image sequence acquisition. According to the absolute phase diagram sequence, the 3D coordinates are calculated to form a 3D image sequence of the human chest and abdomen surface based on the triangulation principle. The sequence expresses the 3D shape of the human chest and abdomen surface at each sampling moment during respiratory movement.

3. Folded Phase Extraction Method

In this paper, 3D-FFA combines 3D Fourier transform with 3D Gauss filtering in the frequency domain to achieve folded phase extraction.

3.1. Folded Phase Extraction Principle

Taking an image in the R fringe image sequence as an example, the principle of folded phase extraction is as follows. Firstly, the intensity of fringe image sequences ir(x,y,t) at different times t can be described as:
i r ( x , y , t ) = a r ( x , y , t ) + b r ( x , y , t ) cos [ 2 π ( f x 0 x + f y 0 y + f t 0 t ) + φ r ( x , y , t ) ]
where x represents the row coordinate of stripe images, y denotes the column coordinate, a r ( x , y , t ) is the background light intensity, b r ( x , y , t ) is the modulation of fringes, f x 0 , f y 0 , and f t 0 are the carrier frequencies in the direction of x , y , and t , and φ r ( x , y , t ) is the phase distribution function. Equation (1) can be further expressed as:
i r ( x , y , t ) = a r ( x , y , t ) + d r ( x , y , t ) exp [ j 2 π ( f x 0 x + f y 0 y + f t 0 t ) ] + d r * ( x , y , t ) exp [ j 2 π ( f x 0 x + f y 0 y + f t 0 t ) ]
where:
d r ( x , y , t ) = 1 2 b r ( x , y , t ) exp [ j φ r ( x , y , t ) ]
d r * ( x , y , t ) = 1 2 b r ( x , y , t ) exp [ j φ r ( x , y , t ) ]
After 3D Fourier transform of Equation (2), we can obtain:
I r ( f x ,   f y ,   f t ) = A r ( f x ,   f y ,   f t ) + D r ( f x f x 0 ,   f y f y 0 ,   f t f t 0 ) + D r * ( f x + f x 0 ,   f y + f y 0 ,   f t + f t 0 )
where f x 0 , f y 0 , and f t 0 are the frequency domain variables in the direction of axes x , y , and t , respectively, A r ( f x , f y , f t ) is the background light spectrum, and D r ( f x f x 0 ,   f y f y 0 ,   f t f t 0 ) and D r * ( f x + f x 0 ,   f y + f y 0 ,   f t + f t 0 ) are the spectra of deformed fringes.
In addition, a 3D filter is used to separate the first level spectrum of D r ( f x f x 0 ,   f y f y 0 ,   f t f t 0 ) and move it to the origin of the frequency domain. After obtaining D r ( f x ,   f y ,   f t ) , the 3D inverse Fourier transform is performed, and the phase distribution function is as follows,
φ r ( x , y , t ) = tan 1 Im { d r ( x , y , t ) } Re { d r ( x , y , t ) }
where Im { d r ( x , y , t ) } denotes the imaginary part of d r ( x , y , t ) and Re { d r ( x , y , t ) } is the real part. Similarly, φ g ( x , y , t ) and φ b ( x , y , t ) can be obtained.

3.2. Three-Dimensional Gauss Filter

When extracting the positive first-order spectrum of fringe image sequences, the 3D filter must also have the function of filtering interference, which is important for 3D Fourier analysis. Because of the interference existing in the environment, the measured object and the measurement system, and the spectrum leakage led by signal truncation, the folding phase error will happen, which will cause 3D measurement errors. Currently, 3D-FFA mainly uses two types of 3D filters [17]: one is the 3D rectangular filter, and the other is the 3D Butterworth filter [18,19]. The former has truncation problems and large leakage errors. For phase unwrapping, the latter has many problems such as a complex algorithm, accumulated error, and unreliability. By comparison, the Gauss filter has the advantages of a small ringing effect and a good effect of eliminating spectrum leakage. For the determined chest and abdomen surface of the measured human, the center frequency and the width of filter in the 3D direction are determined by experiments, which can make it have good adaptability and filter performance. Moreover, previous studies proved that the effect of the 2D Gauss filter is better than the Hanning window and rectangular window [20]. In this paper, the 3D Gauss filter is used as follows,
H ( f x , f y , f t ) = e [ ( f x f 0 x ) 2 2 σ x 2 + ( f y f 0 y ) 2 2 σ y 2 + ( f t f 0 t ) 2 2 σ t 2 ]
where f x 0 , f y 0 , and f t 0 denote the center frequency of the direction of the x , y , and t axis, respectively, and σ x , σ y , and σ t represent the filter widths in the three directions, respectively. Figure 2 gives the schematic diagram of filtering in the frequency domain with 1D-FFA, 2D-FFA, and 3D-FFA. As shown in Figure 2, the cut-off frequency in the x-axis, y-axis, and z-axis directions are f x 1 = f x 0 σ x / 2 and f x 2 = f x 0 + σ x / 2 , f y 1 = f y 0 σ y / 2 and f y 2 = f y 0 + σ y / 2 , and f t 1 = f t 0 σ t / 2 and f t 2 = f t 0 + σ t / 2 , respectively. The cut-off frequency value is determined experimentally for each measured object, so that the measured signal passes through as much as possible, and the interference signal passes as little as possible.
In principle, 3D-FFA filters have stronger anti-interference ability than 1D-FFA and 2D-FFA. The 3D shock interference is taken as an example to explain this theory. The spectrum amplitude of the shock interference δ ( x , y , t ) obtained by 3D Fourier transform is one, and its frequency components cover the entire 3D frequency domain. 1D-FFA can only be filtered along the f x axis in the one-dimensional frequency domain, and its pass band is f x 1 < f x < f x 2 . As shown in Figure 2a, it can only filter out interference signals in the one-dimensional frequency domain. 2D-FFA filters along the f x axis and f y axis in the 2D frequency domain, the pass bands are f x 1 < f x < f x 2 and f y 1 < f y < f y 2 . As shown in Figure 2b, it can filter out interference signals in the 2D frequency domain, which can further significantly weaken the interference signal. 3D-FFA filters along the f x , f y , and f t axis in the 3D frequency domain, and the pass bands are f x 1 < f x < f x 2 , f y 1 < f y < f y 2 , and f t 1 < f t < f t 2 , respectively. As shown in Figure 2c, it can filter out interference signals in the 3D frequency domain, which can weaken the interference signal once again.
Due to adding the relationship between adjacent fringe images in time sequence, we can get the Gaussian filter parameters with individual characteristics by designing and optimizing the filter parameters to reduce the ringing effect and spectrum leakage; this can improve the accuracy of wrapped phase extraction. Moreover, 3D-FFA processes all the images at the same time, so it has high efficiency and is suitable for dynamic measurement.

4. Tri-Frequency Time Phase Unwrapping Method

Firstly, positive integers P r , P g , and P b ( P r < P g < P b ) are used to represent the fringe period of the R, G, and B stripe patterns; then r r , r g , and r b are used to denote the folding phases of the R, G, and B stripe images, and they are the solutions of φ r , φ g , and φ b , respectively in the main value intervals, whose value range is [ 0 , P r ] , [ 0 , P g ] , and [ 0 , P b ] . Therefore, there are the following equations.
φ r × P r / ( 2 × π ) = r r ( mod P r )
φ g × P g / ( 2 × π ) = r g ( mod P g )
φ b × P b / ( 2 × π ) = r b ( mod P b )
In these equations, P r , P g , and P b are the three modulus values, and r r , r g , and r b are the corresponding three remainders.
In addition, N represents the distance between the measured pixel point and the phase origin, which also denotes the absolute phase after the unwrapping processing. Let N r = [ N / P r ] , N g = [ N / P g ] , N b = [ N / P b ] , where [   ] represents a rounding down operation, then N can be written as follows:
N = N r P r + r r = N g P g + r g = N b P b + r b
If P r , P g , and P b have the greatest common divisor P , there are Ω r = P r / P , Ω g = P g / P , Ω b = P b / P , Ω = Ω r × Ω g × Ω b , λ r = Ω / Ω r , λ g = Ω / Ω g , and λ b = Ω / Ω b . Assuming that Ω r , Ω g , and Ω b are mutually prime, then λ r and Ω r are relatively prime, that is the modulus inverse λ ¯ r of λ r exists for Ω r , and there is λ ¯ r λ r 1 ( mod Ω r ) . For the same reason, there are λ ¯ g λ g 1 ( mod Ω g ) and λ ¯ b λ b 1 ( mod Ω b ) .
Let:
q r = [ r r / P ] , q g = [ r g / P ] , q b = [ r b / P ]
Then, there is:
r r = q r P + r c , r g = q g P + r c , r b = q b P + r c
where N r c ( mod P ) and N 0 = [ N / P ] , then:
λ ¯ r λ r q r + λ ¯ g λ g q g + λ ¯ b λ b q b = N 0 ( mod Ω )
so:
N = P N 0 + r c
It is obvious that the folding phases of the R, G, and B fringe images, i.e., the remainders r r , r g , and r b , can be used to unwrap the phase and obtain the absolute phase by Equations (12)–(15).
However, the above discussion does not consider the effect of the residual measurement errors. According to Equation (14), if there are measurement errors Δ r r , Δ r g , and Δ r b in the folding phases r r , r g , and r b , the errors Δ q r , Δ q g , and Δ q b belonging to q r , q g , and q b , respectively will be produced, i.e., | Δ q r | 1 , | Δ q g | 1 , and | Δ q b | 1 . In the meantime, ( λ ¯ r λ r q r + λ ¯ g λ g q g + λ ¯ b λ b q b ) will also produce coarse errors, i.e., | Δ q r λ ¯ r λ r | λ ¯ r × Ω g × Ω b , | Δ q g λ ¯ g λ g | λ ¯ g × Ω r × Ω b , and | Δ q b λ ¯ b λ b | λ ¯ b × Ω r × Ω g . All these can lead the N 0 and N to produce coarse errors Δ N M A X Δ N 0 M A X . Under the condition of λ ¯ r 3 , λ ¯ g 3 , and λ ¯ b 3 , there is | Δ N M A X | 3 × Ω r × Ω g . Consider that the pixel resolution of the existing device is usually limited within Ω b 9 , i.e., | Δ N M A X | Ω / 3 │. Δ N M A X may further cause ( λ ¯ r λ r q r + λ ¯ g λ g q g + λ ¯ b λ b q b ) to make errors in the remainder operation with the modulo Ω = Ω r × Ω g × Ω b , which may result in a larger absolute phase error. To sum up, Δ N M A X can cause a large error in the measurement or even lead the measurement to fail. Therefore, some measures should be adopted to avoid the appearance of Δ N M A X or reduce its influence.
To avoid the appearance of Δ N M A X , assume that the remainder measurement error satisfies Condition A,
| Δ r r / P | < 0.25 , | Δ r g / P | < 0.25 , | Δ r b / P | < 0.25
The measured values r ^ r , r ^ g , and r ^ b of the remainders r r , r g , and r b can be expressed as:
r ^ r / P = [ r ^ r / P ] + { r ^ r / P } , r ^ g / P = [ r ^ g / P ] + { r ^ g / P } , r ^ b / P = [ r ^ b / P ] + { r ^ b / P }
where {   } is a fractional operation. Then, if the difference between the fractional part of the remainder measurement satisfies Condition B, then:
| Δ r r g / P | = | { r ^ r / P } { r ^ g / P } | < 0.5 | Δ r r b / P | = | { r ^ r / P } { r ^ b / P } | < 0.5 | Δ r g b / P | = | { r ^ g / P } { r ^ b / P } | < 0.5
Then, according to the following equation, q r , q g , and q b can be obtained as follows:
q r = [ r ^ r / P ] , q g = [ r ^ g / P ] , q b = [ r ^ b / P ]
Otherwise, according to the following formula, q r , q g , and q b can be calculated:
q r = [ r ^ r / P + 0.5 ] , q g = [ r ^ g / P + 0.5 ] , q b = [ r ^ b / P + 0.5 ]
This ensures that Δ q r = 0 , Δ q g = 0 , Δ q b = 0 , then Δ N 0 = 0 , which can eliminate the coarse error Δ N M A X .
According to Equation (14), N 0 is obtained from q r , q g , and q b , then the absolute phase can be obtained by the following equation:
N = P N 0 + { r ^ r / P } + { r ^ g / P } + { r ^ b / P } 3 P
According to the above equation, the absolute phase error Δ N = Δ r r + Δ r g + Δ r b 3 does not exceed the folding phase errors Δ r r , Δ r g , and Δ r b .
If the residual measurement error does not satisfy Condition A, Δ N M A X may occur. In order to reduce the influence of Δ M A X , the proposed method in this paper increases the part of judging and processing the absolute phase. After the phase unwrapping is completed, the difference Δ N k = N k N k 1 between the absolute phase measurement value N k of each pixel k and the absolute phase measurement value N k 1 of its neighboring pixel is calculated. If it meets the following Condition C,
| Δ N k | < 3 × Ω r × Ω g
We regarded N i as the valid measurement value, and the absolute phase error does not exceed the folding phase errors Δ r r , Δ r g , and Δ r b ; otherwise, | Δ N k | 3 × Ω r × Ω g │, which means the spatial distance between adjacent pixels is no less than 1/3 of the range. This is obviously unreasonable for a relatively flat surface such as the chest and abdomen surface of the human body. If N i is invalid, the pixel can be rejected as an immeasurable point.
If necessary, the interpolation method can be used to obtain the absolute phase based on the surrounding pixels of the absolute phase. Elimination or interpolation usually has little effect on the measurement because the sample points in the image are large and dense. Conversely, if a larger absolute phase error cannot be identified and eliminated, it will affect the measurement result seriously or even make the measurement result unusable. This is also a challenging problem in the 3D measurement of Fourier fringe analysis [21]. To sum up, under Conditions of A, B, and C, Equations (19)–(21) are combined to form our proposed 3D measurement method.
In addition, compared with the three-frequency differential method, the proposed method has a larger unwrapping range. Let P r = P 0 W r , P g = P 0 , P b = P 0 + W b , where W r and W b are positive integers, and the phase unwrapping range of the proposed method is P r g b O = P r P g P b . When phase unwrapping is performed by the three-frequency differential method, the light stripes R and G are used for phase unwrapping to form a synthetic light stripe RG with the phase unwrapping range P r g = P r P g / ( P g P r ) ; then, use the light stripes G and B for phase unwrapping to form a synthetic light stripe GB with a phase unwrapping range P g b = P g P b / ( P b P g ) ; moreover, the phase unwrapping is further performed by using the synthesized light stripes RG and GB. When P g b > P r g , its phase unwrapping range is P r g b H = P r P g P b / ( 2 W r W b + ( W r W b ) P 0 ) . When W r W b , 2 W r W b + ( W r W b ) P 0 2 , then P r g b O 2 P r g b H . Only when W r < W b , 2 W r W b + ( W r W b ) P 0 = 1 , that is P r g b O = P r g b H is possible. According to the pixel resolution of the currently available digital pattern projection device and digital image acquisition device, the phase unwrapping range should be 300 to 10,000 pixels. Only in these two cases, one being P r = 9   pixel , P g = 11   pixel , P b = 14   pixel , P r g b H = 1386   pixels and the other being P r = 12   pixel , P g = 17   pixel , P b = 29   pixel , P r g b H = 5916   pixels , there will be P r g h O = P r g b H . Under these situations, it is difficult to achieve comprehensive optimization of the measurement range, resolution, and anti-interference ability by flexible selection of P r , P g , and P b .

5. Experimental Results and Analysis

5.1. Simulation Experiments of Folding Phase Extraction

Ethical approval to undertake this project was examined by the Human Research Ethics Committee for Non-Clinical Faculties, School of Measurement-Control Technology and Communication Engineering, Harbin University of Science and Technology on 1 March 2019. The title of the project is “Projection on Patient Body Surface in Invasive Surgeries (National Natural Science Foundation of China, 61671190)”. Informed consent form was obtained from the subject. Simulation experiments were conducted by using a planar cosine light image sequence with a size of 768 × 768 pixels and a period P 0 of 35 pixels. To simulate human respiratory movement that approximates periodic motion, let the measured plane do relative paralleled movement to the image plane of the camera 2mm each time and perform a periodic reciprocating motion with a period of 20mm, then collect the plane cosine light images from 120 positions.
i m ( x , y , t ) = 128 + 75 × cos [ 2 π x / P 0 + 2 π ( t 1 ) / 20 ]   m = 1 , 2 , 3 , ,   120 .
The 3D representation of a cosine fringe image sequence with one period in the t-axis direction is shown in Figure 3a. Extract the folding phases of the sixth frame image by 1D-FFA, 2D-FFA, and 3D-FFA, respectively. Then, these extracted folded phases are subtracted from the folded phase, respectively. The folding phase errors of the three methods are shown respectively in Figure 3a–c.
According to Figure 3, in the extracted folded phase by the three methods, there were errors caused by truncation and sampling. The calculated folding phase value of measured points and their standard values were used to calculate the RMS error, and the difference between the maximum folding phase value and the minimum folding phase value was defined as peak-to-valley (PV) error. As shown in Table 1, the PV and RMS of the folded phase error were approximately equal, and they were so small that could be neglected. It can be seen that these three methods could extract the folding phase effectively.
To evaluate the anti-interference ability of 3D-FFA, the following interference signals were added to the cosine fringe image sequence:
γ = 75 × I p × [ 2 × rand ( 768 , 768 , 120 ) 1 ]
where rand() is a function of generating random numbers in [ 0 , 1 ] and I p is the ratio percentage of the interference signal amplitude to the cosine modulation.
By analyzing the fringe images collected in the experiments, the results showed that the noise mainly consisted of salt and pepper noise and Gaussian noise, whose probability distribution curves were superimposed and integrated to form a uniform noise probability curve. Therefore, the uniformly distributed random noise could be used in the simulation experiment, which could fully simulate the effect of noise and was better than the direct superposition of Gaussian noise and salt and pepper noise.
Taking I p = 60 % as an example, the folding phase error of the sixth frame image is shown in Figure 4. It can be seen that the folding phase error extracted by 3D-FFA was obviously smaller than that of the other two methods, which showed that the anti-interference ability of the 3D-FFA was the strongest.
When I p was 1%, 2%, 3%, 5%, 10%, 20%, 40%, 60%, 80%, and 100%, respectively, the folded phase error curves of the three methods are shown in Figure 5. As far as the RMS error and peak-valley error of folded phase were concerned, one was that they increased with the increase of interference; the other was that the results obtained by the 2D-FFA method were significantly lower than the 1D-FFA method; the third was that the results obtained by 3D-FFA method were significantly lower than the 2D-FFA method; the fourth was that the 3D-FFA method and the 2D-FFA method were almost invariant and approximate to the interference when the interference ratio was less than 10%.
In order to compare the anti-interference ability of the three methods quantitatively, the ratio of the folding phase errors of the three methods with different interference percentages is shown in Figure 6. In Figure 6a, the red line is the RMS error ratio of 1D-FFA and 2D-FFA, which was between four and six; the blue line is the RMS error ratio of 2D-FFA and 3D-FFA, which was between one and two. In Figure 6b, the red line is the PV error ratio of 1D-FFA and 2D-FFA, which was between four and eight; the blue line is the PV error ratio of 2D-FFA and 3D-FFA, which was between one and two. Obviously, the anti-interference capability of 2D-FFA method was much better than that of the 1D-FFA method, and the one of the 3D-FFA method was much better than that of the 2D-FFA method; the bigger the interference, the more obvious the superiority of the anti-interference capability. When the interference reached 40%, the anti-interference ability of the 3D-FFA method was about twice that of the 2D-FFA method.

5.2. Chest Model Measurement Experiments

Based on the proposed method in this paper, a 3D experimental apparatus was constructed for human chest and abdomen surface measurement. The device used a projector (InFocus IN82, InFocus Corporation, Wilsonville, OH, USA) to project a color stripe pattern with a resolution of 1024 × 768 pixels. The pattern parameter was P r = 25   pixel , P g = 30   pixel , P b = 35   pixel , and the measured surface stripe images with a resolution of 1624 × 1236 pixels were collected by using a 3CCD industrial camera (AT-200GE, JAI Ltd., Copenhagen, Denmark).
In the measurement experiments, the chest model simulated respiratory movement and reciprocating motion on the guide rail for 20 mm. For each mobile 2 mm, we collected an image as shown in Figure 7. A total of 120 images was collected and sent to the computer to form a sequence of stripe images. Intercept the stripe image sequence with the size of 768 × 768 pixels from the fringe image sequence, and then get the 3D image sequence of the tested area by using the 3D measurement method. Figure 8a–c shows the sixth frame of 3D images formed by taking folded phases using 1D-FFA, 2D-FFA, and 3D-FFA, respectively. All of them could reproduce the 3D surface of the measured area correctly, which verified the 3D measurement method presented in this paper. The visual effect of the three methods was basically the same, because the measurement errors of the three were basically the same when the interference could be ignored in the darkroom.
In order to verify and compare the anti-interference ability of the 3D measurement method in this paper, we added the interference γ of different I p to the stripe sequence of the chest model. Take the sixth frame stripe image as an example, the measurement results are shown in Figure 9. When the interference was 5%, the error of the three measurements was similar.
Beginning with I p = 10 % , the error of 1D-FFA measurement result increased rapidly, and the errors of the other two measurements also increased. The error of measurement based on 2D-FFA was larger than 3D-FFA in the range of 10%–80% for I p , although it was difficult to observe visually. From the beginning of I p = 80 % , the error of measurement based on 2D-FFA was larger than that based on 3D-FFA. It showed that the proposed method had the strongest anti-interference ability, and the greater the interference, the more obvious the advantage.

5.3. Measurement Experiments of the Human Chest and Abdomen Surface

With the proposed 3D measurement method, a 3D measurement experiment of the chest and abdomen surface was conducted during its respiratory movement, as shown in Figure 10a, and the 3D surface depth image sequence of the measured area is shown in Figure 10b. The sixth frame in the 3D image sequence is shown in Figure 10c. That is to say that, the proposed method could reconstruct the 3D surface of the human chest and abdomen correctly. In Figure 10c, the enlarged part of the square area on the left is shown on the right, which indicated that there was maximum measurement error in the indoor light environment, which formed an undetectable blank spot.
In the measurement process, interference was added artificially. When the interference was small, the measurement results were basically unchanged. When the interference reached 30%, the measurement results are shown in Figure 10d, which indicated that the undetectable points increased significantly, and the measurement results started to become significantly worse. If necessary, the measured value of the white area could be obtained by the difference between the effective measurement values of its peripheral adjacent points.
However, for a relatively flat surface such as the surface of the human chest and abdomen, eliminating the isolated white spot or obtaining the blank area by the difference had little effect on the measurement results. As shown in Figure 11, a fixed point near the diaphragm was selected, and the height (z) of the point in the image sequence changed with time (t) to form the height (z)-time (t) respiratory movement curve of the point. The measurement results of two typical respiratory movement states, deep breathing and rapid breathing, are given in Figure 11a,b, respectively.
The measured results showed that the proposed method in this paper could recover the 3D shape of the surface of the chest and abdomen at different times correctly, realize the dynamic 3D measurement of the surface, and obtain the respiratory motion trajectory of the surface points. This method had strong anti-interference ability and could eliminate the maximum measurement error or its influence basically.

6. Conclusions

In this paper, a 3D dynamic surface measurement method was proposed for the human chest and abdomen. The RGB trichromatic cosine stripe pattern was synthesized into a color stripe projection pattern. The projection pattern was invariable, and one measurement could be realized by collecting one stripe image. It had less pattern projection and image acquisition time, which provided the foundation for dynamic measurement. Color correction and color separation for the color stripe image sequence were used to form the RGB monochromatic deformation stripe image sequence with a different period. The 3D Fourier transform and the 3D Gaussian filter were combined to carry out 3D-FFA to extract the folding phase. The interference effect could be reduced, and the measurement accuracy could be improved by increasing the time dimension frequency domain filter. The phase unwrapping operation depended on the difference between the fractional parts of the folded phase error, which could ensure that the absolute phase error did not exceed the measurement error of the remainder. Removal or interpolation based on absolute phase measurement could eliminate or reduce the effect of gross absolute phase error. The phase unwrapping based on the remainder equation had a larger phase unwrapping range and preferred space.
In view of this proposed method, the theoretical analysis and experimental results showed that the method presented in this paper could obtain the 3D image sequence of the human chest and abdomen surface and the respiratory motion curve of the chest and abdomen surface points in the respiratory movement effectively and accurately. It had the characteristics of strong anti-interference ability and a wide range of development and could eliminate gross absolute phase error. It should be pointed out that although the method in this paper achieved the analysis of five frames of 3D images per second and respiratory motion in about four seconds, its operation time should be reduced from both the perspective of improving the analysis effect and engineering demand. Therefore, we plan to adopt multithread optimization measures to improve the operation speed.

Author Contributions

This article was completed by all authors. H.W. designed and implemented the classification algorithm. S.Y. made the experimental analysis of the algorithm. X.Y. participated in the writing of the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (61671190, 61571168), the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (UNPYSCT-2017086), and Fundamental Research Foundation for Universities of Heilongjiang Province (LGYC2018JQ014).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fassi, A.; Schaerer, J.; Fernandes, M.; Riboldi, M.; Sarrut, D.; Baroni, G. Tumor tracking method based on a deformable 4D CT breathing motion model driven by an external surface surrogate. Int. J. Radiat. Oncol. Biol. Phys. 2014, 88, 182–188. [Google Scholar] [CrossRef] [PubMed]
  2. Fayad, H.; Pan, T.; Pradier, O.; Visvikis, D. Patient specific respiratory motion modeling using a 3D patient’s external surface. Med. Phys. 2012, 39, 3386–3395. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Povsic, K.; Jezersek, M.; Mozina, J. Real-time 3D visualization of the thoraco-abdominal surface during breathing with body movement and deformation extraction. Physiol. Meas. 2015, 36, 1497–1516. [Google Scholar] [CrossRef] [PubMed]
  4. Xiao, H.; Hou, L.; You, Y. A Measurement Method of Multi-line Structure light Stereo vision based on Image Fusion. J. Sichuan Univ. (Eng. Sci. Ed.) 2015, 47, 154–158. [Google Scholar]
  5. Sam, V.D.J.; Dirckx, J.J.J. Real-time structured light profilometry: A review. Opt. Lasers Eng. 2016, 87, 18–31. [Google Scholar]
  6. Zuo, C.; Tao, T.; Feng, S.; Huang, L.; Asundi, A.; Chen, Q. Micro Fourier Transform Profilometry (mu FTP): 3D shape measurement at 10,000 frames per second. Opt. Lasers Eng. 2018, 102, 70–91. [Google Scholar] [CrossRef] [Green Version]
  7. Lu, L.; Ding, Y.; Luan, Y.; Luan, Y.; Yin, Y.; Liu, Q.; Xi, J. Automated approach for the surface profile measurement of moving objects based on PSP. Opt. Express 2017, 25, 32120–32131. [Google Scholar] [CrossRef] [Green Version]
  8. Zhong, M.; Su, X.Y.; Chen, W.J.; You, Z.S.; Lu, M.T.; Jing, H.L. Modulation measuring profilometry with auto-synchronous phase shifting and vertical scanning. Opt. Express 2014, 22, 31620–31634. [Google Scholar] [CrossRef]
  9. EI-Haddad, M.T.; Tao, Y.K.K. Automated stereo vision instrument tracking for intraoperative OCT guided anterior segment ophthalmic surgical maneuvers. Biomed. Opt. Express 2015, 6, 3014–3031. [Google Scholar] [CrossRef] [Green Version]
  10. Su, X.Y.; Zhang, Q.C. Dynamic 3-D shape measurement method: A review. Opt. Lasers Eng. 2010, 48, 191–204. [Google Scholar] [CrossRef]
  11. Shirato, H.; Shimizu, S.; Kitamura, K.; Nishioka, T. Four-dimensional treatment planning and fluoroscopic real-time tumor tracking radiotherapy for moving tumor. Int. J. Radiat. Oncol. Biol. Phys. 2000, 48, 435–442. [Google Scholar] [CrossRef]
  12. Cao, S.P.; Cao, Y.P.; Zhang, Q.C. Fourier transform profilometry of a single-field fringe for dynamic objects using an interlaced scanning camera. Opt. Commun. 2016, 367, 130–136. [Google Scholar] [CrossRef]
  13. Geng, J. Structured-light 3D surface imaging: A tutorial. Adv. Opt. Photonics 2011, 3, 128–160. [Google Scholar] [CrossRef]
  14. Padilla, M.; Servin, M.; Garnica, G. Fourier analysis of RGB fringe-projection profilometry and robust phase-demodulation methods against crosstalk distortion. Opt. Express 2016, 24, 15417–15428. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Caspi, D.; Kiryati, N.; Shamir, J. Range imaging with adaptive color structured light. IEEE Trans. Pattern Anal. Mach. Intell. 1998, 20, 470–480. [Google Scholar] [CrossRef]
  16. Hu, Y.; Chen, Q.; Zhang, Y.Z.; Feng, S.J.; Tao, T.Y.; Li, H.; Yin, W.; Zuo, C. Dynamic microscopic 3D shape measurement based on marker-embedded Fourier transform profilometry. Appl. Opt. 2018, 57, 772–780. [Google Scholar] [CrossRef]
  17. Shi, H.J.; Zhu, F.P.; He, X.Y. Low-Frequency Vibration Measurement Based on Spatiotemporal Analysis of Shadow Moire. Acta Opt. Sin. 2011, 31, 120–124. [Google Scholar]
  18. Abdul-Rahman, H.S.; Gdeisat, M.A.; Burton, D.R.; Lalor, M.J.; Lilley, F.; Abid, A. Three-dimensional Fourier fringe analysis. Opt. Lasers Eng. 2008, 46, 446–455. [Google Scholar] [CrossRef]
  19. Zhang, Q.C.; Hou, Z.L.; Su, X.Y. 3D fringe analysis and phase calculation for the dynamic 3D measurement. AIP Conf. Proc. 2010, 1236, 395–400. [Google Scholar]
  20. Bu, P.; Chen, W.J.; Su, X.Y. Analysis on measuring accuracy of Fourier transform profilometry due to different filtering window. Laser J.. 2003, 24, 43–45. [Google Scholar]
  21. Zheng, D.L.; Da, F.P.; Kemao, Q.; Seah, H.S. Phase-shifting profilometry combined with Gray-code patterns projection: Unwrapping error removal by an adaptive median filter. Opt. Express 2017, 25, 4700–4713. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Schematic diagram of the 3D measurement principle and system of the human chest and abdomen surface. (1) Pattern projection and image acquisition. (2) Image color correction and separation. (3) Folded phase extraction. (4) Folded phase unwrapping. (5) Three-dimensional image sequence acquisition.
Figure 1. Schematic diagram of the 3D measurement principle and system of the human chest and abdomen surface. (1) Pattern projection and image acquisition. (2) Image color correction and separation. (3) Folded phase extraction. (4) Folded phase unwrapping. (5) Three-dimensional image sequence acquisition.
Sensors 20 01091 g001
Figure 2. Schematic diagram of filtering in the frequency domain with three methods. FFA, Fourier fringe analysis.
Figure 2. Schematic diagram of filtering in the frequency domain with three methods. FFA, Fourier fringe analysis.
Sensors 20 01091 g002
Figure 3. Cosine light image sequence and folding image error of one image.
Figure 3. Cosine light image sequence and folding image error of one image.
Sensors 20 01091 g003
Figure 4. Folding phase error of a stripe image after adding interference.
Figure 4. Folding phase error of a stripe image after adding interference.
Sensors 20 01091 g004
Figure 5. The folded phase error curves of the sixth frame stripe images with different interference percentages. (a) Comparison of the root mean square of the folded phase error; (b) comparison of peak and valley values of folded phase errors.
Figure 5. The folded phase error curves of the sixth frame stripe images with different interference percentages. (a) Comparison of the root mean square of the folded phase error; (b) comparison of peak and valley values of folded phase errors.
Sensors 20 01091 g005
Figure 6. The ratio curves between the folded phase errors of the three methods with different interference ratios of the sixth frame. (a) RMS ratio curve of folded phase error; (b) PV ratio curve of folded phase error.
Figure 6. The ratio curves between the folded phase errors of the three methods with different interference ratios of the sixth frame. (a) RMS ratio curve of folded phase error; (b) PV ratio curve of folded phase error.
Sensors 20 01091 g006
Figure 7. Image of the chest model after the projected stripe pattern.
Figure 7. Image of the chest model after the projected stripe pattern.
Sensors 20 01091 g007
Figure 8. Measurement results of the sixth frame stripe image of the chest model.
Figure 8. Measurement results of the sixth frame stripe image of the chest model.
Sensors 20 01091 g008
Figure 9. Measurement results of the sixth stripe image of the chest model with different interference.
Figure 9. Measurement results of the sixth stripe image of the chest model with different interference.
Sensors 20 01091 g009aSensors 20 01091 g009b
Figure 10. Surface of the human chest and abdomen and its measurement results. (a) Measured surface of the chest and abdomen; (b) 3D surface depth image sequence of the measured area; (c) measurement results of the sixth stripe image and the local magnification map when ignoring the interference; (d) measurement results of the sixth stripe image and the local magnification map after adding 30% interference.
Figure 10. Surface of the human chest and abdomen and its measurement results. (a) Measured surface of the chest and abdomen; (b) 3D surface depth image sequence of the measured area; (c) measurement results of the sixth stripe image and the local magnification map when ignoring the interference; (d) measurement results of the sixth stripe image and the local magnification map after adding 30% interference.
Sensors 20 01091 g010aSensors 20 01091 g010b
Figure 11. The height (z)-time (t) respiratory motion curve at a point near the diaphragm. (a) Curve of the point’s motion in deep breathing; (b) curve of the point’s motion in rapid breathing.
Figure 11. The height (z)-time (t) respiratory motion curve at a point near the diaphragm. (a) Curve of the point’s motion in deep breathing; (b) curve of the point’s motion in rapid breathing.
Sensors 20 01091 g011
Table 1. Folding phase error of the stripe image (unit: rad).
Table 1. Folding phase error of the stripe image (unit: rad).
Methods1D-FFA2D-FFA3D-FFA
PV0.01560.01560.0145
RMS0.00400.00400.0044

Share and Cite

MDPI and ACS Style

Wu, H.; Yu, S.; Yu, X. 3D Measurement of Human Chest and Abdomen Surface Based on 3D Fourier Transform and Time Phase Unwrapping. Sensors 2020, 20, 1091. https://doi.org/10.3390/s20041091

AMA Style

Wu H, Yu S, Yu X. 3D Measurement of Human Chest and Abdomen Surface Based on 3D Fourier Transform and Time Phase Unwrapping. Sensors. 2020; 20(4):1091. https://doi.org/10.3390/s20041091

Chicago/Turabian Style

Wu, Haibin, Shuang Yu, and Xiaoyang Yu. 2020. "3D Measurement of Human Chest and Abdomen Surface Based on 3D Fourier Transform and Time Phase Unwrapping" Sensors 20, no. 4: 1091. https://doi.org/10.3390/s20041091

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop