CN111275669A - Priori information guided four-dimensional cone beam CT image reconstruction algorithm - Google Patents

Priori information guided four-dimensional cone beam CT image reconstruction algorithm Download PDF

Info

Publication number
CN111275669A
CN111275669A CN202010032744.7A CN202010032744A CN111275669A CN 111275669 A CN111275669 A CN 111275669A CN 202010032744 A CN202010032744 A CN 202010032744A CN 111275669 A CN111275669 A CN 111275669A
Authority
CN
China
Prior art keywords
image
phase
motion
prior
reconstruction
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010032744.7A
Other languages
Chinese (zh)
Other versions
CN111275669B (en
Inventor
牟轩沁
职少华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202010032744.7A priority Critical patent/CN111275669B/en
Publication of CN111275669A publication Critical patent/CN111275669A/en
Application granted granted Critical
Publication of CN111275669B publication Critical patent/CN111275669B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

The invention discloses a four-dimensional cone-beam CT image reconstruction algorithm guided by prior information, the prior image reconstructed by the invention does not need to introduce other imaging modality images by utilizing the existing data, and the image quality is high and is easy to obtain; the method can produce high-quality phase resolution images without depending on other imaging modes such as clinical CT and the like. The linear estimation model provided by the invention is simple and effective, and can accurately estimate the motion compensation time sequence image, accurately depict the motion track of respiration, and simultaneously retain the detail information of the image. The method can solve various 4D CBCT cases with large breathing movement, and a large number of experimental results verify the effectiveness of the method; according to the invention, from the novel point of view of importance of the initial image to the final result in the four-dimensional cone beam CT iterative reconstruction framework, the motion compensation time sequence image is constructed to have good space-time resolution, and the iterative reconstruction image can be guided to gradually approach to the expected image quality.

Description

Priori information guided four-dimensional cone beam CT image reconstruction algorithm
Technical Field
The invention belongs to the technical field of X-ray imaging, and particularly relates to a four-dimensional cone-beam CT image reconstruction algorithm guided by prior information.
Background
Radiation therapy is one of the main approaches to the treatment of malignant tumors (commonly known as cancer). Cone Beam CT (CBCT) is widely used in the field of radiotherapy, particularly in Image Guided Radiation Therapy (IGRT). The cone beam CT has the main functions of positioning a patient, positioning a treatment target according to a reconstruction result and providing real-time update of patient information for a radiotherapy process. Generally, the detector of cone-beam CT in this field has high spatial resolution, which is very convenient for head positioning, but faces a great challenge for lung imaging. In the cone beam CT imaging time of 1-4 minutes per week, in the relatively long scanning process, respiration, heartbeat and irregular gastrointestinal movement of a patient inevitably occur, and a lung tissue structure with serious movement blurring can be reconstructed, so that the accurate positioning of a target area in radiotherapy is influenced. From the above, one problem with cone-beam CT that has not been fully solved to date is the problem of imaging dynamic objects.
Among the problems of cone-beam CT imaging for dynamic targets, the time-Dimension cone-beam CT technique (4D CBCT) which incorporates time factors into reconstruction was first proposed in 2005 and is also called Four-dimensional cone-beam CT. The main strategy of the technology is to synchronously acquire corresponding respiratory motion signals in the cone beam CT scanning process, perform time phase classification on projection images according to the respiratory motion signals, respectively reconstruct images under the current phase, and finally obtain a group of time sequence cone beam CT images. This set of images may provide motion trajectory information for the time dimension of the moving organ, with high temporal resolution. Although motion artifacts in the reconstruction results are significantly suppressed at different breathing phases, the projection data assigned at each phase is equivalent to sparsely sampling the projection data for all angles. In general, the motion breathing phase in one cycle can be time-divided into ten states, and the number of projections corresponding to each phase state is one tenth of all the acquired projection data, and this sampling form cannot satisfy Shannon-Nyquist sampling theorem. On the other hand, the projection data in each phase are in a beam-bunching distribution, which further aggravates the serious sparse angle problem. Therefore, undersampled streak artifacts are necessarily introduced in the reconstructed image, while the spatial resolution of the reconstructed image is severely reduced.
In recent years, a great deal of research work has emerged aiming at improving the image quality of four-dimensional cone-beam CT, and reconstructing high-quality time series images while accurately estimating and locating a dynamic target tumor in a respiratory state. Current mainstream solutions can be divided into: the image quality is improved by estimating a motion change model between different phases, and four-dimensional CT iterative reconstruction is carried out by utilizing prior information such as image correlation between phases. The first method considers that images with different phases can be non-rigidly registered by calculating a Deformable Vector Field (DVF) in combination with a prior image to obtain images with different phases. However, this method has a high requirement on the accuracy of the registration algorithm and is prone to introduce new error points. In order to improve the registration accuracy, some methods need to introduce high-quality CT images of other modalities as a reference. The second strategy is a phase-dependent four-dimensional cone-beam CT reconstruction algorithm, which utilizes the high correlation of the patient at different respiratory phases to perform spatial-Temporal joint constraint to eliminate the bar artifacts in the reconstructed image, such as Temporal Non-local mean filtering (TNLM) and Temporal-spatial motion compensation type total variation. However, the capability of such methods for image artifact removal is still increasing.
Disclosure of Invention
The invention aims to overcome the defects and provide a four-dimensional cone beam CT image reconstruction algorithm guided by prior information, which can generate high-quality phase resolution images.
In order to achieve the above object, the present invention comprises the steps of:
acquiring original projection data P in cone beam CT imaging equipment, extracting breathing signals of cone beam CT projection images, and performing phase distribution on projection data to obtain a projection data subset Pt；
Step two, utilizing complete projection data P and projection data subset PtCarrying out reconstruction to obtain a prior image XpriorAnd a set of three-dimensional volume images containing T phases, a priori image XpriorAnd the four-dimensional cone beam CT image sequence is Xt，t＝1,…,T；
Respectively extracting background information and motion characteristics of each phase reconstruction image according to the correlation of the phase image sequence;
fourthly, the motion characteristics of each phase reconstruction image and the prior image XpriorCombining, building linear estimation model, and constructing motion compensation image under each phase
Figure BDA0002364919510000031
Step five, the motion compensation images under each phase position
Figure BDA0002364919510000032
As the initial value of the iterative reconstruction method, the final image is reconstructed
Figure BDA0002364919510000033
The objective function is:
Figure BDA0002364919510000034
wherein A istRepresenting the slave image X at the t-th phasetTo the corresponding phase projection data PtIs a system matrix of |)TVA regularization constraint term representing the total variation, and λ represents a regularization factor.
In the first step, when respiratory signals are extracted, the number of respiratory states in a respiratory cycle is set to be T, and the number of respiratory phases is obtained.
If the respiratory signals cannot be synchronously extracted, collecting by an Amsterdam smooth method, wherein the Amsterdam smooth method comprises the following specific steps:
firstly, performing point-by-point derivation on a projection image obtained by a detector at each projection angle along the projection angle, and summing and averaging pixels in a row along the horizontal direction of the derived image to obtain a column of images;
secondly, arranging the obtained column images according to a projection angle to generate an AS image;
and thirdly, carrying out derivation on the obtained AS image along the horizontal direction, and extracting a respiratory signal from the derived AS image.
In the first step, projection data is carried outWhen the phase is distributed, firstly, all the original projection data P are distributed to the corresponding breathing state according to the phase information of the breathing signal to obtain T groups of projection data, and the projection data under any phase T is Pt。
In the second step, an objective function is established:
Figure BDA0002364919510000041
respectively reconstructing prior images X according to different projection measurement data corresponding to the target functionpriorAnd four-dimensional cone-beam CT image sequence Xt；
In the objective function, a represents a system matrix from an image X to be reconstructed to corresponding projection data P, r (X) represents a regularization term of the reconstructed image, and λ represents a regularization factor.
In step three, the four-dimensional cone beam CT image sequence X obtained by reconstruction is repeatedtExpressing, extracting a low-rank component and a sparse component in each phase image, and constructing a model function as follows:
Figure BDA0002364919510000042
wherein L istRepresenting the t-th phase image XtLow rank component of (2), StRepresenting the t-th phase image XtThe motion component of (a); r is a regularization parameter, | Lt‖*A representation matrix LtRepresents the sum of its singular values, | St‖1The 1 norm of the matrix is represented and s.t. is the mathematical notation abbreviation of the constraint.
In the fourth step, the constructed linear estimation model is a motion blur artifact estimation model, and a high-quality initialization image in the iterative reconstruction process is constructed.
The concrete method of the step four is as follows:
first, calculating the motion characteristics S of each phasetRemoving fuzzy artifacts in the prior image to obtain a synthetic background image XsynBackground image XsynThe method does not contain motion components under each phase, has clear background detail organization structure, and has the calculation formula as follows:
Figure BDA0002364919510000043
second, the motion information of the current phase is compensated to the synthesized background image XsynConstructing a group of motion compensation image sequences with clear background and reflecting motion characteristics under each phase
Figure BDA0002364919510000044
Figure BDA0002364919510000051
Compensating the motion of the image
Figure BDA0002364919510000052
And (4) as an initial input value of the minimization target function, circularly iterating the solved image through the SART algorithm.
Compared with the prior art, the reconstructed prior image does not need to introduce other imaging modality images by utilizing the existing data, and has high image quality and easy acquisition; the method can produce high-quality phase resolution images without depending on other imaging modes such as clinical CT and the like. The linear estimation model provided by the invention is simple and effective, and can accurately estimate the motion compensation time sequence image, accurately depict the motion track of respiration, and simultaneously retain the detail information of the image. The method can solve various 4D CBCT cases with large breathing movement, and a large number of experimental results verify the effectiveness of the method; the method is based on the novel angle that the importance of the initial image to the final result in the four-dimensional cone-beam CT iterative reconstruction framework is taken as a basis, the motion compensation time sequence image is constructed to have good space-time resolution, the iterative reconstruction image can be guided to gradually approach to the expected image quality, and the method has important guiding and inspiring significance for the future four-dimensional cone-beam CT reconstruction.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a general schematic of phase assignment of projection data;
FIG. 3 is a graphical illustration of the number of projection data assignments versus relative window width and step size;
FIG. 4 is a diagram of an intermediate result display for an XCAT die body embodiment of the present invention;
FIG. 5 is a graph of experimental results comparing different methods in an XCAT die embodiment; wherein, (a) is a cross section and coronal section slice image of a true value image and a reconstructed prior image under a phase [ 0% 10% ] and a phase [ 40% 50% ]; (b) cross-sectional versus coronal slice results plots of the results reconstructed at phase [ 0% 10% ] for five different methods; (c) cross-sectional versus coronal slice result plots of the results reconstructed at phase [ 40% 50% ] for six different methods;
FIG. 6 is a cross-sectional and coronal slice result plot of the reconstruction results of different methods at phase [ 0% 10% ];
FIG. 7 is a cross-sectional and coronal slice result diagram of five different phase reconstruction results obtained by the method of the present invention.
Detailed Description
The invention will be further explained with reference to the drawings.
Referring to fig. 1, the present invention comprises the steps of:
acquiring original projection data P in cone beam CT imaging equipment, extracting breathing signals of cone beam CT projection images, and performing phase distribution on projection data to obtain a projection data subset Pt(ii) a When respiratory signals are extracted, the number of respiratory states in a respiratory cycle is set to be T, and the number of respiratory phases is obtained; if the respiratory signals cannot be synchronously extracted, acquiring by an Amsterdam method, wherein the Amsterdam method comprises the following specific steps:
firstly, performing point-by-point derivation on a projection image obtained by a detector at each projection angle along the projection angle, and summing and averaging pixels in a row along the horizontal direction of the derived image to obtain a column of images;
secondly, arranging the obtained column images according to a projection angle to generate an AS image;
and thirdly, in order to highlight the motion signal effect, derivation is carried out on the obtained AS image along the horizontal direction, and the respiratory signal is extracted from the derived AS image.
When the phase distribution of the projection data is carried out, firstly, all the original projection data P are distributed to the corresponding breathing state according to the phase information of the breathing signal to obtain T groups of projection data, and the projection data under any phase T is Pt. Where phase information refers to partitioning into T groups of equal time intervals for each cycle in the motion signal, each interval being considered a relatively "stationary" state of respiratory motion.
Step two, utilizing complete projection data P and projection data subset PtCarrying out reconstruction to obtain a prior image XpriorAnd a set of three-dimensional volume images containing T phases, a priori image XpriorAnd the four-dimensional cone beam CT image sequence is XtT is 1, …, T; establishing an objective function:
Figure BDA0002364919510000061
respectively reconstructing prior images X according to different projection measurement data corresponding to the target functionpriorAnd four-dimensional cone-beam CT image sequence Xt；
In the objective function, a represents a system matrix from an image X to be reconstructed to corresponding projection data P, r (X) represents a regularization term of the reconstructed image, and λ represents a regularization factor. In order to optimize the minimization function, the method uses a joint iterative Reconstruction algorithm (SART) to solve. For the regularization term, the method adopts a regularization method of Total Variation Minimization constraint (TV) to inhibit and smooth noise and bar artifacts in an image to be reconstructed, namely R (X) | X |TVThe calculation formula is as follows:
Figure BDA0002364919510000071
in the formula, i, j and k represent pixel coordinates of the image in three directions of x, y and z;
respectively extracting background information and motion characteristics of each phase reconstruction image according to the correlation of the phase image sequence; for the four-dimensional cone beam CT image sequence X obtained by reconstructiontExpressing, extracting a low-rank component and a sparse component in each phase image, and constructing a model function as follows:
Figure BDA0002364919510000072
wherein L istRepresenting the t-th phase image XtLow rank component of (2), StRepresenting the t-th phase image XtThe motion component of (a); r is a regularization parameter, | Lt‖*A representation matrix LtRepresents the sum of its singular values, | St‖1Representing the 1 norm of the matrix. This model is a Robust Principal Component Analysis model (RPCA). The invention solves the RPCA model by using a PCP (Principal Component Pursuit) method, finally reduces the dimension of a four-dimensional cone beam CT image sequence, and reconstructs an image X under each phasetDecomposed into a highly correlated background matrix LtAnd a sparse matrix S containing motion feature informationt(ii) a The sparse matrix StIncluding motion state information at the current phase and a small amount of artifacts and noise, StReferred to as motion components.
Fourthly, the motion characteristics of each phase reconstruction image and the prior image XpriorCombining, building linear estimation model, and constructing motion compensation image under each phase
Figure BDA0002364919510000073
The constructed linear estimation model is a motion blurring artifact estimation model, and a high-quality initialization image in the iterative reconstruction process is constructed.
First, calculating the motion characteristics S of each phasetRemoving fuzzy artifacts in the prior image to obtain a synthetic background image XsynBackground image XsynThe method does not contain motion components under each phase, has clear background detail organization structure, and has the calculation formula as follows:
Figure BDA0002364919510000081
second, the motion information of the current phase is compensated to the synthesized background image XsynConstructing a group of motion compensation image sequences with clear background and reflecting motion characteristics under each phase
Figure BDA0002364919510000082
Figure BDA0002364919510000083
Compensating the motion of the image
Figure BDA0002364919510000084
And (3) as an initial input value of a minimization target function, circularly iterating through an SART algorithm to obtain a solved image, wherein the minimization target function is as follows:
Figure BDA0002364919510000085
step five, the motion compensation images under each phase position
Figure BDA0002364919510000086
As the initial value of the iterative reconstruction method, the final image is reconstructed
Figure BDA0002364919510000087
The objective function is:
Figure BDA0002364919510000088
wherein A istRepresenting the slave image X at the t-th phasetTo the corresponding phase projection data PtIs a system matrix of |)TVA regularization constraint term representing the total variation, and λ represents a regularization factor.
Example (b):
referring to fig. 1,2 and 3, the present invention comprises the steps of:
step one, projection data phase classification;
and obtaining original complete projection data P of the CT imaging equipment at an angle of 0-360 degrees. In the embodiment, the real-time respiratory motion signal is obtained in the imaging process, and the respiratory signal extraction is not needed. Next, the respiratory motion signal is divided into T phases at equal time intervals, where T is 10 in this embodiment. As shown in fig. 2, corresponding projections are extracted and assigned to the same group according to the same breathing phase time t, and finally 10 sets of projection data sets P at different phases are obtainedtAnd t is 1,2, …, 10. In order to avoid the too small number of projection angles in each phase, the invention introduces the phase width and the step length, and the projection data in the adjacent phases are distributed and partially overlapped, thereby ensuring that the number of the projection angles in each phase is increased while the time resolution of the four-dimensional cone beam CT is not reduced.
Reconstructing different phase images and prior images;
constructing an iterative reconstruction framework based on TV minimization:
Figure BDA0002364919510000091
specifically, an SART method is adopted for updating, and a TV-based minimization regular term is carried out to suppress noise in a reconstructed image. The TV calculation formula is:
Figure BDA0002364919510000092
in the formula, i, j, k represents pixel coordinates of the image in three directions of x, y, z.
1) Reconstructing a prior image;
iteratively reconstructing a prior image X based on the reconstruction framework using the complete projection data Pprior。
Figure BDA0002364919510000093
Because of more projection angles, in order to accelerate reconstruction, the projection data of each angle is divided into a plurality of subsets (Ordered subsets, OS) at equal intervals, and each Subset data is updated sequentially every iteration. Meanwhile, the non-negative constraint of the reconstructed image is considered in each step of iterative updating process. X in FIG. 4priorCross-sectional slice and coronal slice results representing a priori reconstructed images reconstructed by the above-described method.
2) Reconstructing different phase images;
the projection data P after phase classificationt( t 1,2, …,10), and the reconstructed images X at different phases are independently reconstructedt(T ═ 1, …, T); for the current phase t, its objective function becomes:
Figure BDA0002364919510000094
wherein A istRepresenting images X from which images are to be reconstructedtTo corresponding projection data PtA system matrix under a geometric projection relationship; at the same time, the physical meaning of the linear attenuation coefficient of a substance is considered, which is constantly positive, so that each pixel in the reconstructed image is non-negatively constrained in an iterative process. PRI (Phase-resolved Image, PRI) in FIG. 4 represents [ 40% -50% of the total amount of the solution reconstructed by the above method]Cross-sectional slices and coronal slices at phase.
3) RPCA-based motion component extraction;
aiming at a group of three-dimensional image sequences obtained by reconstruction in the step 1), wherein X is ═ X1,X2,…,Xt,…,XT]，Xt∈RI×J×KT is 1, …, T. Using different phasesThe high correlation between the two matrixes converts the similar structure into a subspace matrix L, so that the subspace matrix L has a relatively low rank. Meanwhile, the different motion state information S ═ S reflected by each phase image is separated1,S2,…,St,…,ST](StMotion information representing the t-th phase). The final overall optimization objective function is:
min‖L‖*+r‖S‖1s.t.X＝L+S
wherein | L |*Representing the kernel norm of the matrix L, representing the sum of its singular values, and calculated as | L |*＝∑j＝1σj(Li)；‖S‖1Representing the 1 norm of the matrix by the formula
Figure BDA0002364919510000101
r > 0 is a constant. This model is a Robust Principal Component Analysis model (RPCA).
Due to the huge number of voxels of the reconstructed three-dimensional image, the three-dimensional image sequence needs to be processed independently layer by layer in the actual calculation process. Is provided with
Figure BDA0002364919510000102
(N ═ I × J) is the kth layer (0) in the tth phase<K is less than or equal to K) column vectorization results of the two-dimensional image, and the results of the same kth layer in T different phases are combined into a matrix Mk∈RN×TThen, the RPCA modeling is performed, i.e. the new sub-targeting function is:
Figure BDA0002364919510000103
wherein L iskRepresents MkLow rank component of (2), SkThe motion information component of the kth layer in the tth phase image; at this time, the parameter takes the value of
Figure BDA0002364919510000104
In the invention, R is subjected to principal component tracking (PCP) methodCarrying out optimization solution on the PCA model, and finally carrying out reconstruction on an image X under each phasetDecomposed into a highly correlated background matrix LtAnd a sparse matrix S containing motion feature informationt. The sparse matrix StMainly contains motion state information at the current phase and a small amount of artifacts and noise. The cross-sectional slice and coronal slice results of the background matrix L and the motion features obtained by modeling calculation by the above method are shown in fig. 4, respectively.
4) Constructing a motion compensation initial image sequence;
the prior image X obtained by the calculation of the step 2) and the step 3) is usedpriorWith the motion characteristics S at different phasest(T-1, …, T) are combined to construct a high quality initialization image in an iterative reconstruction process based on a proposed motion blur artifact estimation model.
First, motion component information S in each phase is calculatedtTo remove the fuzzy artifact in the prior image and obtain a synthetic background image Xsyn. The background image is characterized by not containing motion components under each phase and having a clear background detail organization structure. The calculation formula is as follows:
Figure BDA0002364919510000111
x in FIG. 4synThe cross-section slice and coronal plane slice results of the synthesized background image calculated by the motion blur artifact estimation model are shown.
Second, compensating the motion component information of the current phase to the synthesized background image to construct a group of motion compensation image sequence with clear background and reflecting the motion characteristics of each phase
Figure BDA0002364919510000112
Figure BDA0002364919510000113
The motion compensated image in fig. 4 represents the cross-sectional slice versus coronal slice results at [ 40% -50% ] phase, as compared to the true-valued image.
5) Individual phase image update reconstruction (motion compensated image guided);
compensating the motion of the image
Figure BDA0002364919510000114
The solved image is iterated through the SART algorithm loop as an initial input value based on the TV minimization objective function. The minimization objective function is:
Figure BDA0002364919510000115
the iterative reconstruction update steps based on TV minimization for step 2) and step 5) are as follows:
initializing a reconstructed image:
Figure BDA0002364919510000121
setting parameters: the number of phases T is 10;
total number of cycles N main10; number of fidelity item subcycles Niter＝5；
Iterative gradient descent step size factor λ1＝1；
Number of TV regularization term loops NTV＝5；
Gradient descent step factor lambda in TV2＝0.002；
Figure BDA0002364919510000122
The effectiveness of the prior information guided four-dimensional cone-beam CT image reconstruction method is verified through a large number of simulation data experiments and actual data experiments by combining specific experiment results.
The following experiments were specifically performed:
a) preparing data: a set of phantom data and a set of actual patient projection data. Wherein, XCAT software is used for simulating and generating the die bodies containing the detailed anatomical structure, and each die body contains 10 motion states. Specifically, the phantom 1 has a small breathing amplitude during one cycle, but has more detailed information (such as blood vessels, abdominal soft tissues, etc.). For the actual data, we acquired projection data of the chest region. The cone-beam CT used is a flat panel imaging device integrated on a TrueBeam medical accelerator (Varian medical systems, Palo Alto, CA). The measured geometric parameters and reconstruction parameters of the simulated projections and actual data of the phantom are shown in table 1.
TABLE 1 reconstruction parameters of XCAT phantom and actual thoracic data
Figure BDA0002364919510000131
b) The superiority of the proposed algorithm is verified by comparing different four-dimensional cone beam CT reconstruction methods; wherein the comparison method comprises PRI, PICCS, RPCA-4DCT and SMART method;
c) detailed quantitative analysis is carried out on simulation data by using indexes such as Root Mean Square Error (RMSE), Entropy function (Entropy), Normalized Mutual Information value (NMI) and the like;
specifically, the root mean square error is used to measure the absolute error between the result reconstructed by the adopted method and the true image, and the lower the RMSE value is, the smaller the absolute error is. The calculation formula is as follows:
Figure BDA0002364919510000132
in the formula, the first step is that,
Figure BDA0002364919510000133
representing the results obtained by different methods of reconstruction,
Figure BDA0002364919510000134
representing the corresponding true value image; n is a radical ofi,j,kRepresenting the number of all voxels of the three-dimensional image;
the entropy function mainly measures the degree of the adopted method for restraining the bar-shaped artifacts, and the larger the entropy value is, the lower the degree of the adopted method for restraining the image artifacts is; the normalized mutual information represents the correlation between the reconstructed image X and the true value image Y, and the higher the NMI value is, the higher the reconstruction quality is. The calculation formula is as follows:
Figure BDA0002364919510000141
Figure BDA0002364919510000142
wherein, X represents a reconstructed image, GT represents a true value image (GT); h (q) a gray histogram representing different pixel intensities, N representing the sum of the number of pixels; h (GT) is the entropy of the truth image GT, H (X, GT) represents the joint entropy of the reconstructed image X and the truth image GT, and NMI (X, GT) is the normalized mutual information value of the reconstructed image X and the truth image GT.
The simulation data reconstruction result is shown in fig. 5, the actual patient data reconstruction result is shown in fig. 6 and 7, compared with other methods, the method can reconstruct high-quality images with different phases, effectively inhibits the strip artifacts and the noise caused by the sparse angle problem, and meanwhile, the method can recover the detail structure information of the sectional images, such as tiny blood vessels (fig. 5) and focuses (fig. 6 and 7) of the lung position.
In terms of quantitative analysis, table 2 gives the root mean square error values of XCAT reconstruction results for different methods at [ 0% 10% ] and [ 40% 50% ] phases. Table 3 shows the entropy function value and mutual information value of the XCAT simulation reconstruction result by different methods. As can be seen from the table, compared to other methods, the method of the present invention reconstructs an image having the smallest root mean square error value, the lowest entropy function value as shown in table 3, and the highest mutual information value as shown in table 3.
TABLE 2 RMS error values of XCAT reconstruction results at [ 0% 10% ] and [ 40% 50% ] phases for different methods
PRI SMART RPCA-4DCT PICCS Motion compensated image Results of the proposed method
[40％50％]Phase position 148.01 134.77 137.38 130.31 124.82 93.73
[0％-10％]Phase position 151.52 117.98 137.94 127.19 115.37 86.90
TABLE 3 entropy function value and mutual information value of XCAT simulation reconstruction result by different methods
Figure BDA0002364919510000143
Figure BDA0002364919510000151
The method obviously inhibits the bar artifact caused by the sparse angle problem of the four-dimensional cone beam CT while eliminating the motion artifact. The method explains the problem of four-dimensional cone-beam CT reconstruction from the novel point of view of constructing a high-quality initialization image to guide the reconstruction process. Firstly, reconstructing a prior image by utilizing all collected measured projection data; then, extracting the motion characteristics of each phase image by using a principal component analysis model; then, combining the motion characteristics with the prior image to construct a linear estimation model to obtain a group of high-quality motion compensation time sequence images; the motion compensation image is used as an initial image of a four-dimensional cone-beam CT iterative reconstruction frame, and finally a group of four-dimensional cone-beam CT images with high space-time resolution are obtained. The invention has the following characteristics: 1) the prior image utilizes the existing data without introducing other imaging modality images, so that the image quality is high and the image is easy to obtain; 2) the proposed linear estimation model is simple and effective, and the estimated motion compensation time sequence image can accurately depict the motion track of respiration, and meanwhile, the detail information of each phase image is retained.

Claims (8)

1. A priori information guided four-dimensional cone-beam CT image reconstruction algorithm is characterized by comprising the following steps:
acquiring original projection data P in cone beam CT imaging equipment, extracting breathing signals of cone beam CT projection images, and performing phase distribution on projection data to obtain a projection data subset Pt；
Step two, utilizing complete projection data P and projection data subset PtCarrying out reconstruction to obtain a prior image XpriorAnd a set of three-dimensional volume images containing T phases, a priori image XpriorAnd the three-dimensional volume image form a four-dimensional cone beam CT, and the four-dimensional cone beam CT image sequence isXt，t＝1,…,T；
Respectively extracting background information and motion characteristics of each phase reconstruction image according to the correlation of the phase image sequence;
fourthly, the motion characteristics of each phase reconstruction image and the prior image XpriorCombining, building linear estimation model, and constructing motion compensation image under each phase
Figure FDA0002364919500000012
Step five, the motion compensation images under each phase position
Figure FDA0002364919500000013
As the initial value of the iterative reconstruction method, the final image is reconstructed
Figure FDA0002364919500000014
The objective function is:
Figure FDA0002364919500000011
wherein A istRepresenting the slave image X at the t-th phasetTo the corresponding phase projection data PtIs a system matrix of |)TVA regularization constraint term representing a total variation, λ represents a regularization factor,
Figure FDA0002364919500000015
is the initial value of the image to be reconstructed.
2. The a priori information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 1, wherein in the first step, when the respiratory signal is extracted, the number of respiratory states in a respiratory cycle is set to T, that is, the number of respiratory phases is obtained.
3. The prior information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 2, wherein if the respiratory signals cannot be synchronously extracted, the acquisition is performed by an Amsterdam method, and the Amsterdam method comprises the following specific steps:
firstly, performing point-by-point derivation on a projection image obtained by a detector at each projection angle along the projection angle, and summing and averaging pixels in a row along the horizontal direction of the derived image to obtain a column of images;
secondly, arranging the obtained column images according to a projection angle to generate an AS image;
and thirdly, carrying out derivation on the obtained AS image along the horizontal direction, and extracting a respiratory signal from the derived AS image.
4. The prior information guided four-dimensional cone-beam CT image reconstruction algorithm according to claim 1, wherein in the first step, when performing phase assignment of projection data, all original projection data P are first assigned to corresponding respiratory states according to phase information of respiratory signals to obtain T sets of projection data, and projection data at any phase T is Pt。
5. The a priori information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 1, wherein in step two, an objective function is established:
Figure FDA0002364919500000021
respectively reconstructing prior images X according to different projection measurement data corresponding to the target functionpriorAnd four-dimensional cone-beam CT image sequence Xt；
In the objective function, a represents a system matrix from an image X to be reconstructed to corresponding projection data P, r (X) represents a regularization term of the reconstructed image, and λ represents a regularization factor.
6. The a priori information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 1, wherein the step of reconstructing the four-dimensional cone-beam CT image is performed in stepIn the third step, the four-dimensional cone beam CT image sequence X obtained by reconstruction is reconstructedtExpressing, extracting a low-rank component and a sparse component in each phase image, and constructing a model function as follows:
Figure FDA0002364919500000022
wherein L istRepresenting the t-th phase image XtLow rank component of (2), StRepresenting the t-th phase image XtThe motion component of (a); r is a regularization parameter, | Lt‖*A representation matrix LtRepresents the sum of its singular values, | St‖1Representing the 1 norm of the matrix.
7. The a priori information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 1, wherein in step four, the constructed linear estimation model is a motion blur artifact estimation model, and a high quality initialization image in an iterative reconstruction process is constructed.
8. The a priori information guided four-dimensional cone-beam CT image reconstruction algorithm of claim 7, wherein the specific method of step four is as follows:
first, calculating the motion characteristics S of each phasetRemoving fuzzy artifacts in the prior image to obtain a synthetic background image XsynBackground image XsynThe method does not contain motion components under each phase, has clear background detail organization structure, and has the calculation formula as follows:
Figure FDA0002364919500000031
second, the motion information of the current phase is compensated to the synthesized background image XsynConstructing a group of motion compensation image sequences with clear background and reflecting motion characteristics under each phase
Figure FDA0002364919500000033
Figure FDA0002364919500000032
Compensating the motion of the image
Figure FDA0002364919500000034
And (4) as an initial input value of the minimization target function, circularly iterating the solved image through the SART algorithm.
CN202010032744.7A 2020-01-13 2020-01-13 Priori information guided four-dimensional cone beam CT image reconstruction algorithm Active CN111275669B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010032744.7A CN111275669B (en) 2020-01-13 2020-01-13 Priori information guided four-dimensional cone beam CT image reconstruction algorithm

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010032744.7A CN111275669B (en) 2020-01-13 2020-01-13 Priori information guided four-dimensional cone beam CT image reconstruction algorithm

Publications (2)

Publication Number Publication Date
CN111275669A true CN111275669A (en) 2020-06-12
CN111275669B CN111275669B (en) 2022-04-22

Family

ID=71000184

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010032744.7A Active CN111275669B (en) 2020-01-13 2020-01-13 Priori information guided four-dimensional cone beam CT image reconstruction algorithm

Country Status (1)

Country Link
CN (1) CN111275669B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111951346A (en) * 2020-08-18 2020-11-17 安徽工程大学 4D-CBCT reconstruction method combining motion estimation and space-time tensor enhancement expression
CN112819911A (en) * 2021-01-23 2021-05-18 西安交通大学 Four-dimensional cone beam CT reconstruction image enhancement algorithm based on N-net and CycN-net network structures
CN113189543A (en) * 2021-04-27 2021-07-30 哈尔滨工程大学 Interference suppression method based on motion compensation robust principal component analysis
CN113450345A (en) * 2021-07-19 2021-09-28 西门子数字医疗科技(上海)有限公司 Image processing method, image processing device, electronic equipment and storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130002659A1 (en) * 2010-02-12 2013-01-03 The Regents Of The University Of California Graphics processing unit-based fast cone beam computed tomography reconstruction
CN104321773A (en) * 2012-03-31 2015-01-28 瓦里安医疗***公司 4D cone beam CT using deformable registration
CN104899907A (en) * 2015-07-06 2015-09-09 东南大学 Sparse angle CT image reconstruction method based on gamma prior
US20160174921A1 (en) * 2014-12-19 2016-06-23 Ion Beam Applications S.A. Method and imaging system for determining a reference radiograph for a later use in radiation therapy
CN109035360A (en) * 2018-07-31 2018-12-18 四川大学华西医院 A kind of compressed sensing based CBCT image rebuilding method
CN110148215A (en) * 2019-05-22 2019-08-20 哈尔滨工业大学 A kind of four-dimensional MR image reconstruction method based on smoothness constraint and local low-rank restricted model
CN110390361A (en) * 2019-07-25 2019-10-29 安徽工程大学 A kind of 4D-CBCT imaging method based on motion compensation study

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130002659A1 (en) * 2010-02-12 2013-01-03 The Regents Of The University Of California Graphics processing unit-based fast cone beam computed tomography reconstruction
CN104321773A (en) * 2012-03-31 2015-01-28 瓦里安医疗***公司 4D cone beam CT using deformable registration
US20160174921A1 (en) * 2014-12-19 2016-06-23 Ion Beam Applications S.A. Method and imaging system for determining a reference radiograph for a later use in radiation therapy
CN104899907A (en) * 2015-07-06 2015-09-09 东南大学 Sparse angle CT image reconstruction method based on gamma prior
CN109035360A (en) * 2018-07-31 2018-12-18 四川大学华西医院 A kind of compressed sensing based CBCT image rebuilding method
CN110148215A (en) * 2019-05-22 2019-08-20 哈尔滨工业大学 A kind of four-dimensional MR image reconstruction method based on smoothness constraint and local low-rank restricted model
CN110390361A (en) * 2019-07-25 2019-10-29 安徽工程大学 A kind of 4D-CBCT imaging method based on motion compensation study

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
HAO YAN ET AL: "A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging", 《MEDICAL PHYSICS》 *
LUDWIG RITSCHL ET AL: "Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatio-temporal sparsity prior", 《PHYSICS IN MEDICINE AND BIOLOGY》 *
SHAOHUA ZHI ET AL: "Artifacts reduction method for phase-resolved Cone-Beam CT (CBCT) images via a prior-guided CNN", 《SPIE MEDICAL IMAGING》 *
徐霄: "基于形变配准的4D-CBCT稀疏角度重建研究", 《中国优秀硕士学位论文全文数据库（电子期刊）医药卫生科技辑》 *
杨轩等: "基于运动补偿的压缩感知4D-CBCT优质重建", 《J SOUTH MED UNIV》 *
陈美玲等: "运动配准先验图像的4D-CBCT优质重建", 《J SOUTH MED UNIV》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111951346A (en) * 2020-08-18 2020-11-17 安徽工程大学 4D-CBCT reconstruction method combining motion estimation and space-time tensor enhancement expression
CN112819911A (en) * 2021-01-23 2021-05-18 西安交通大学 Four-dimensional cone beam CT reconstruction image enhancement algorithm based on N-net and CycN-net network structures
CN112819911B (en) * 2021-01-23 2022-10-25 西安交通大学 Four-dimensional cone beam CT reconstruction image enhancement algorithm based on N-net and CycN-net network structures
CN113189543A (en) * 2021-04-27 2021-07-30 哈尔滨工程大学 Interference suppression method based on motion compensation robust principal component analysis
CN113450345A (en) * 2021-07-19 2021-09-28 西门子数字医疗科技(上海)有限公司 Image processing method, image processing device, electronic equipment and storage medium

Also Published As

Publication number Publication date
CN111275669B (en) 2022-04-22

Similar Documents

Publication Publication Date Title
CN111275669B (en) Priori information guided four-dimensional cone beam CT image reconstruction algorithm
Cai et al. Cine cone beam CT reconstruction using low-rank matrix factorization: algorithm and a proof-of-principle study
Terpstra et al. Deep learning-based image reconstruction and motion estimation from undersampled radial k-space for real-time MRI-guided radiotherapy
Isola et al. Fully automatic nonrigid registration‐based local motion estimation for motion‐corrected iterative cardiac CT reconstruction
CN106846430B (en) Image reconstruction method
Hinkle et al. 4D CT image reconstruction with diffeomorphic motion model
CN110390361B (en) 4D-CBCT imaging method based on motion compensation learning
CN115605915A (en) Image reconstruction system and method
Chen et al. A new variational model for joint image reconstruction and motion estimation in spatiotemporal imaging
CN110570489B (en) Motion compensation high-quality 4D-CBCT image reconstruction method based on bilateral filtering design
US10388036B2 (en) Common-mask guided image reconstruction for enhanced four-dimensional cone-beam computed tomography
Jung et al. Deep learning cross-phase style transfer for motion artifact correction in coronary computed tomography angiography
Karakatsanis et al. Quantitative PET image reconstruction employing nested expectation-maximization deconvolution for motion compensation
CN112819911B (en) Four-dimensional cone beam CT reconstruction image enhancement algorithm based on N-net and CycN-net network structures
Huang et al. A biomechanical modeling-guided simultaneous motion estimation and image reconstruction technique (SMEIR-Bio) for 4D-CBCT reconstruction
CN105976412B (en) A kind of CT image rebuilding methods of the low tube current intensity scan based on the sparse regularization of offline dictionary
Liu et al. Low-dose CBCT reconstruction via 3D dictionary learning
Gupta et al. Neural Computed Tomography
Hinkle et al. 4D MAP image reconstruction incorporating organ motion
CN109767410B (en) Lung CT and MRI image fusion algorithm
CN111951346B (en) 4D-CBCT reconstruction method combining motion estimation and space-time tensor enhancement representation
Liu et al. 4D-CBCT reconstruction via motion compensataion learning induced sparse tensor constraint
CN115049752A (en) PET respiratory motion image artifact registration correction method based on three-dimensional convolutional neural network
Dhou et al. Motion-based projection generation for 4D-CT reconstruction
Wang et al. An unsupervised dual contrastive learning framework for scatter correction in cone-beam CT image

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant