Adaptive direction and step size multi-spectral ct iterative reconstruction with scatter correction

Abstract

The application discloses a kind of self-adapting direction and step of multi-energy spectrum CT iterative reconstruction method with scattering correction, comprising: step 1, setting multi-energy spectrum CT projection reconstruction model with scattering correction;Step 2, using iterative algorithm to solve out the increment of base material distribution projection value angle by angle, update the current base material image to new base material image;Wherein, through self-adapting direction and step to quickly obtain convergent base material projection image.The scheme provided by the application is used to remove scattering, and the result of base material decomposition is more accurate, and the reconstruction result artifact is less, which lays a foundation for quantitative analysis.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of computed tomography (CT), and in particular to a multi-spectral CT iterative reconstruction method with scatter correction and adaptive direction and step length. BACKGROUND

[0002] Computed tomography (CT) is a technique that reconstructs the internal structure of an object without destroying the object by collecting a series of X-ray attenuation data at different angles, and then reconstructing the internal structure of the object by a specific algorithm and a computer. CT imaging has many advantages, such as non-contact, no image overlap, high resolution, easy processing, storage and transmission, and has been widely used in many fields such as medicine, industry, biology, petroleum, materials, security inspection, etc.

[0003] Traditional CT image reconstruction models and algorithms are based on the assumption of single-energy rays, and the reconstructed image is the spatial distribution of the linear attenuation coefficient of the measured object at a certain energy. However, the X-ray generated by the X-ray machine used in the CT widely used in clinical and industrial fields is multi-energy, and the corresponding CT image is actually reconstructed by a set of scanning data of rays at a certain energy spectrum through a single-energy CT image reconstruction algorithm. The reconstructed CT image can be approximately regarded as the spatial distribution of the linear attenuation coefficient of the measured object corresponding to X-rays at a certain energy. Due to the inconsistency between the image reconstruction model and the physical process of data acquisition, combined with factors such as photon scattering, detector efficiency, and signal-to-noise ratio of data, the CT value of the reconstructed image will be distorted in some way, resulting in homoplasy (the same material corresponding to different CT values) and heteroplasy (different materials corresponding to the same CT value). Typical image CT value distortion includes "hardening artifact", "ring artifact", "photon starvation artifact", "scattering artifact", "half-volume effect artifact", etc. Therefore, the material differentiation ability of conventional clinical CT and industrial CT is relatively limited.

[0004] Dual-energy CT has a significantly higher material differentiation ability than traditional CT, and has achieved a series of successful applications in medical and industrial CT imaging. In dual-energy CT imaging, the linear attenuation coefficient of the object is approximately expressed as a linear combination of two basis functions, and by obtaining two sets of X-ray scanning data corresponding to two energy spectra, the coefficient distribution of the two basis functions (i.e. basis images) is reconstructed. Through the two basis images, the equivalent atomic number and equivalent electron density of the measured object can be obtained, and thus the material can be distinguished and quantitative analysis can be performed.

[0005] The dual-energy CT image reconstruction problem can be attributed to the solution of a class of nonlinear integral equations in mathematics, that is, two basis images are reconstructed from two sets of nonlinear projections corresponding to two energy spectrums. There are two types of approximate representations of the linear attenuation coefficient of the measured object. One type is to select two known materials and only rely on the linear attenuation coefficient of X-ray energy as the basis function (commonly known as basis material decomposition). The other type is to select the attenuation coefficient of X-ray photoelectric absorption and the attenuation coefficient of Compton scattering as the basis function (commonly known as effect decomposition). The dual-energy CT image reconstruction methods can be divided into three categories: pre-processing method, post-processing method and direct iterative reconstruction method. Taking the basis material decomposition as an example to illustrate the differences between the three methods, the effect decomposition is similar. The pre-processing method is to first decompose the material in the projection domain to obtain the projection data of the corresponding basis material image, and then reconstruct the basis material image according to the conventional CT reconstruction algorithm. This method usually uses known energy spectrum to perform forward projection on different thicknesses of basis materials, establishes a lookup table of the relationship between thickness and projection, or uses a high-order polynomial to establish the relationship between the line integral of the attenuation projection and the basis image, and then uses different thickness combinations of the basis material to calibrate and fit. This method is only suitable for the case where the projection data of different energy spectrums is compatible, such as using a double-layer detector or a photon counting detector to collect energy spectrum projection data. The post-processing method refers to first reconstructing the nonlinear projection data under each energy spectrum respectively, and then realizing material decomposition in the image domain. The classic representative of this method is the EDEC empirical method, the MDIR method, and various algorithms developed based on photon counting detectors. The direct iterative reconstruction method directly solves the basis image from the nonlinear integral equation through iteration. The representative of this method is ASD-POCS, ASD-NC-POCS and the iterative solution method based on linearization approximation of the nonlinear integral equation. The typical representative of the iterative solution method based on linearization approximation of the nonlinear integral equation is E-ART algorithm, E-SART algorithm, IFBP algorithm and OPMT which is a tilt correction algorithm extended and accelerated from E-ART.

[0006] However, the existing approximate iterative solution method based on linearization of the nonlinear integral equation has the problems of slow convergence rate, weak noise resistance, inadaptation to inconsistent ray paths, and ignoring the influence of scattering signals on multi-energy CT images. SUMMARY

[0007] The purpose of the present application is to provide a multi-energy CT iterative reconstruction method with scatter correction and adaptive direction and step length to overcome or at least alleviate at least one of the above-mentioned defects of the prior art.

[0008] To achieve the above-mentioned purpose, the present application provides a multi-energy CT iterative reconstruction method with scatter correction and adaptive direction and step length, comprising:

[0009] Step 1, the multi-energy CT projection reconstruction model with scattering correction is expressed by the following formula:

[0010]

[0011] in: For the multicolor projection data of X-rays with scattering correction along path l at the k-th energy spectrum, F(l) = (F1(l), F2(l), ..., F M (l)) τ =(∫ l f1(x)dl,∫ l f2(x)dl,…,∫ l f M (x)dl) τ , representing the projection values ​​of the distribution of M base materials along the ray path l, F m (l)=∫ l f m (x)dl, m=1,2,…,M, M is the number of base materials, f m (x) represents the density distribution of the m-th base material in the measured object, and τ represents the transpose of the vector; E min and E max These are the lowest and highest energies of a photon, respectively; S k (E) is the normalized equivalent energy spectrum of the k-th energy spectrum; μ m (E), m=1,2,…,M is the linear attenuation coefficient of the m-th base material; For the k-th energy spectrum S k The set of ray paths of (E), where K is the number of energy spectra;

[0012] Step 2: Using an iterative algorithm, calculate the increment of the base material distribution projection value angle by angle, and update the current base material image to a new base material image; including:

[0013] Let the current base material distribution projection value be... This represents the estimated projection of the base material distribution after n iterations, where n represents the number of iterations, and is initialized to 0.

[0014] For all rays at the current angle Calculated and scattering intensity Sc l and by order Calculate ΔF (n) (l);

[0015] Let F (n+1) (l)=F (n) (l)+ΔF (n) (l), updated to obtain new base material distribution projection value F (n+1)(l) and back-project to obtain the base material image updated in the n+1th iteration.

[0016] Preferably, ΔF (n) (l) includes:

[0017] The function P k (F (n+1) (l)) is first-order Taylor expanded at the point F (n) (l) to obtain the following expression:

[0018]

[0019] wherein, is the function P k The partial derivative of F(l) at the point F (n) (l) represents the weighted average of the attenuation coefficients of the respective base materials under the current path;

[0020]

[0021] The function P k,F (F (n) (l)) is normalized to obtain

[0022]

[0023] wherein,

[0024] wherein, is the directed distance from the current point F (n) (l) to the hyperplane represented by equation (6);

[0025] ΔF (n) (l) is solved according to equations (6) and (7).

[0026] Preferably, ΔF (n) (l) includes:

[0027] Suppose that there are r compatible measurement data of the ray under the current angle, r≤K, and that equation (7) is used to obtain:

[0028]

[0029]

[0030] wherein equation (9) is the matrix form of equation (8),

[0031] The direction of the minimum solution of equation (9) is taken as the main direction, and solving the minimum solution of equation (9) is equivalent to solving the following equation (10):

[0032]

[0033] where ΔF1 (n) (l) represents that formula (10) is established and the minimum ΔF (n) (l) value;

[0034] The main direction is represented as:

[0035]

[0036] In the auxiliary direction, the following formula (15) and (16) are set:

[0037]

[0038]

[0039] wherein formula (16) is a matrix form of formula (15), Let:

[0040] Take the auxiliary direction as

[0041]

[0042] According to the main direction and the auxiliary direction , an adaptive descending direction

[0043]

[0044] wherein a is a constant greater than or equal to 0;

[0045] Set the length of ΔF (n) (l) is constrained by the farthest distance from the current point to each hyperplane, that is k=1,2,…,K, c is a constant;

[0046] Let the step size s (n) (l) be the weighted average of the absolute distance of F (n) (l) to r hyperplanes, which is expressed by the following formula:

[0047]

[0048] wherein, represents the directed distance of F (n) (l) to the kth hyperplane, and ω k represents the weight of ;

[0049] Adaptive descent direction and step size (n) (l) is obtained

[0050] Delta F (n) (l) = s (n) (l) D (n) (l) …… (20).

[0051] The present application has the following advantages due to the above technical solutions:

[0052] With the scheme provided by the present application, after removing scattering, the decomposition result of the base material is more accurate, the reconstruction result has fewer artifacts, and a foundation is laid for quantitative analysis. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 A flowchart of a self-adaptive direction and step size multi-spectral CT iterative reconstruction method with scatter correction provided by an embodiment of the present application.

[0054] Figure 2 A result diagram of the self-adaptive direction and step size multi-spectral CT iterative reconstruction method with scatter correction provided by an embodiment of the present application applied to a water-based image.

[0055] Figure 3 A result diagram of the self-adaptive direction and step size multi-spectral CT iterative reconstruction method with scatter correction provided by an embodiment of the present application applied to a bone-based image.

[0056] Figure 4 A result diagram of the self-adaptive direction and step size multi-spectral CT iterative reconstruction method with scatter correction provided by an embodiment of the present application applied to a pseudo-single-energy image. DETAILED DESCRIPTION

[0057] In the drawings, the same or similar notations are used to indicate the same or similar elements or elements having the same or similar functions. The embodiments of the present application will be described in detail below with reference to the accompanying drawings.

[0058] In the description of the present application, the terms "center", "longitudinal", "transverse", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", and the like indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present application and simplifying the description, and do not indicate or imply that the devices or elements referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as limiting the scope of protection of the present application.

[0059] The technical features in the embodiments and the implementation forms of the present application can be combined with each other in the case of no conflict, and are not limited to the embodiments or the implementation forms in which the technical features are located.

[0060] The present application will be further described below in conjunction with the drawings and specific embodiments. It should be noted that the technical solutions and design principles of the present application are described in detail below only with one optimized technical solution, but the protection scope of the present application is not limited to this.

[0061] This paper relates to the following terms, in order to facilitate understanding, the meaning is explained as follows. Those skilled in the art should understand that the following terms can also have other names, but any other name should be considered consistent with the terms listed herein without deviating from its meaning.

[0062] The embodiment of the present application provides a multi-spectrum CT iterative reconstruction method with scatter correction of adaptive direction and step length, as shown in Figure 1 , comprising:

[0063] Step 1, the multi-spectrum CT projection reconstruction model with scatter correction is represented by the following formula:

[0064]

[0065] Wherein: is the polychromatic projection data of the X-ray with scatter correction along the path l under the kth energy spectrum, F(l) = (F1(l), F2(l), …, F M (l)) τ = (∫ l f1(x)dl,∫ l f2(x)dl,…,∫ l f M (x)dl) τ , represents the distribution projection value of M base materials along the ray path l, F m (l) = ∫ l f m (x)dl, m = 1, 2, …, M, M is the number of base materials, f m (x) represents the density distribution of the mth base material in the measured object, τ represents the transpose of the vector; E min and E max are the lowest energy and the highest energy of the photons respectively; S k (E) is the normalized equivalent energy spectrum of the kth energy spectrum; μ m (E), m = 1, 2, …, M is the linear attenuation coefficient of the mth base material; is the set of ray paths of the kth energy spectrum S k (E); K is the number of energy spectra.

[0066] where the scatter corrected multi-energy projection model expression (1) is obtained by the following way:

[0067] When the multi-energy X-rays pass through the object, its linear attenuation coefficient is not only related to the measured object, but also related to the energy of the X-rays. Therefore, the X-ray intensity passing through the measured object can be expressed as

[0068]

[0069] where I0(l) represents the initial flow intensity of the X-rays along the path l; I(l) represents the flow intensity after passing through the measured object along the X-ray path l; μ(x, E) is the linear attenuation coefficient of the object at the spatial position x for the X-rays with energy E; S k (E) is the normalized equivalent energy spectrum (related to the ray source, filter, detector crystal, etc.); E min and E max are the lowest energy and the highest energy of the photons, respectively; l is the ray path; Sc l is the scatter signal intensity received by the detector pixel corresponding to the path l; is the set of ray paths for the kth energy spectrum S k (E), and K is the number of energy spectra.

[0070] Considering the influence of the scatter signal on the multi-energy CT image quality, the measured attenuation data is first scatter corrected. After each update of the basis material projection image, the scatter signal intensity Sc l received by the detector pixel corresponding to the ray path l is simulated by the method gQMCFFD based on quasi-Monte Carlo and forced detection combination. Denote the scatter corrected attenuation data as

[0071]

[0072] Taking the logarithmic transformation of the above formula, the following scatter corrected multi-energy projection model

[0073]

[0074] where, is the scatter corrected multi-color projection data along the path l. Based on the basis material decomposition method, μ(x, E) is decomposed into the following form

[0075]

[0076] where μ m (E), m = 1, 2, …, M is the linear attenuation coefficient of the mth basis material, and f m(x), m = 1, 2, …, M is the density distribution of the m-th basis material in the object under test, and M is the number of basis materials. Let F(l) = (F1(l), F2(l), …, F M (l)) τ = (∫ l f1(x)dl,∫ l f2(x)dl,…,∫ l f M (x)dl) τ be the distribution projection values of the M basis materials along the ray path l. Therefore, the multi-energy spectral projection model with scatter correction can be expressed as the following M-element nonlinear equation group

[0077]

[0078] Step 2: Use an iterative algorithm to solve the increment of the basis material distribution projection value for each angle, and update the current basis material image to a new basis material image.

[0079] Before the basis material image f m (x), m = 1, 2, …, M is solved by the adaptive direction and step length multi-energy spectral CT iterative reconstruction method with scatter correction (ADS-SC) proposed in the application, the linear attenuation coefficient μ m (E) of the basis material is given in advance, and after the basis material image f m (x) is updated in each iteration, the scatter signal intensity received by the detector pixel corresponding to the ray path l is simulated, and then the energy spectrum S k (E) is estimated by the EM algorithm with scatter correction. For example, the values of the linear attenuation coefficient μ m (E), m = 1, 2, …, M of the basis material are given in advance by looking up a table, and after the basis material image f m (x), m = 1, 2, …, M is updated in each iteration, the scatter signal intensity Sc l received by the detector pixel corresponding to the ray path l is quickly simulated by using a method combining quasi-Monte Carlo and forced detection (gQMCFFD), and the energy spectrum S k (E), k = 1, 2, …, K is estimated by the EM algorithm with scatter correction, and the detector response function is known.

[0080] Let the current basis material distribution projection value be , which represents the estimated value of the basis material distribution projection after n iterations, and n represents the number of iterations, and n = 0 at initialization.

[0081] For all rays under the current angle, the basis material distribution projection value and the scatter intensity Sc l are calculated, and by letting Calculate ΔF (n) (l) ;

[0082] Further, let F (n+1) (l) = F (n) (l) + ΔF (n) (l), update the new base material distribution projection value F (n+1) (l), and use the SART algorithm to back project to obtain the base material image updated after the n+1 iteration

[0083] Wherein, the process of calculating ΔF (n) (l) includes:

[0084] The first order Taylor expansion is performed on the function P k (F (n+1) (l) ) at the point F (n) (l), to obtain the following expression:

[0085]

[0086] Wherein, is the partial derivative of the function P k about F (l) at F (n) (l), also known as the normal direction of the hyperplane (6), indicating the weighted average value of the attenuation coefficient of each base material under the current path;

[0087]

[0088] The normalization processing is performed on P k,F (F (n) (l) ) to obtain

[0089] Wherein,

[0090] According to equations (6) and (7), ΔF (n) (l) can be solved.

[0091] In equation (7), the geometric meaning of D (n) (l) is the directed distance from the current point F (n) (l) to the hyperplane (the hyperplane represented by equation (6) ).

[0092] Solving ΔF (n) (l) according to equations (6) and (7) includes:

[0093] By adapting the direction D (n) (l) and the step size s (n) (l) of ΔF (n)(l) to quickly obtain the converged substrate projection image F(l)=(F1(l),F2(l),…,F M (l)) τ .

[0094] In this invention, through the main direction and auxiliary directions A combination to achieve ΔF (n) (l) Direction D (n) The adaptation of (l) is given below. The main direction is given below. auxiliary direction and step size s (n) (l) Calculation method.

[0095] Suppose that there are r compatible measurement data for the ray at the current angle, r≤K, then from equation (7) we get:

[0096]

[0097]

[0098] Equation (9) is the matrix form of equation (4).

[0099] Taking the direction of the minimum solution of equation (9) as the principal direction, solving the minimum solution of equation (9) is equivalent to solving the optimization problem shown in equation (10):

[0100]

[0101] Wherein, ΔF1 (n) (l) indicates that equation (10) holds and Minimum ΔF (n) (l) value.

[0102] The main direction is then represented as:

[0103]

[0104] The process of obtaining the expression for the principal direction according to equation (10) includes:

[0105] When matrix When the rank is r, this optimization problem is a quadratic convex optimization with a unique solution. It can be solved by constructing a Lagrangian function:

[0106]

[0107] Where λ=(λ1,λ2,…λ r ) τ By taking the partial derivatives of the function with respect to each variable and setting the derivatives to zero, we can obtain...

[0108]

[0109]

[0110] Equation (12) gives ΔF (n) (l) is substituted into equation (13) to obtain the parameter equations

[0111]

[0112] where is the covariance matrix of is a symmetric positive definite matrix. Solving the equations (14) and substituting into (12) gives the solution of the optimization problem (10) and the principal direction .

[0113] The solution of the auxiliary direction is introduced as follows. At the current point F (n) (l), the direction of each equation is known, i.e. the partial derivative k of the function P (n) (F (n) (l)) is calculable. The difference is that the first r functions with scatter-corrected projection data are known, while the last K-r measurements are unknown. A reasonable idea is to keep the unknown function values unchanged (only move along the tangent plane) when calculating the increment.

[0114] Based on the above analysis, in the present application, the following equations (15) and (16) are set in the auxiliary direction:

[0115]

[0116]

[0117] where equation (16) is the matrix form of equation (15), Let:

[0118] then the auxiliary direction is

[0119]

[0120] where equation (16) can be directly solved, or indefinite, or overdetermined. In the present application, the general solution form of the unconstrained optimization with parameter λ is considered

[0121]

[0122] If the above optimization model is not solved directly, the case of overdetermined equation (16) (m>K) can be solved using the least squares method, or mK equations can be discarded and then solved according to the main direction solution method.

[0123] Then, according to the main direction and auxiliary directions Obtain the adaptive descent direction

[0124]

[0125] Where a is a constant greater than or equal to 0.

[0126] In this invention, the step size s (n) (l) is F (n) (l) The weighted average of the absolute distances to the r hyperplanes is expressed by the following formula:

[0127]

[0128] in, F represents (n) (l) the directed distance to the k-th hyperplane, ω k express The weights;

[0129] Using adaptive descent direction and step size s (n) (l) obtained

[0130] ΔF (n) (l)=s (n) (l)D (n) (l)……(20).

[0131] In this invention, the adaptive step size refers to the correction amount ΔF. (n) The length of (l) can be constrained by the farthest distance from the current point to each hyperplane, i.e. k = 1, 2, ..., K, where c is a constant. In this paper, a step size s is chosen. (n) (l) is the weighted average of the absolute distances to each hyperplane, thus obtaining the above equation (19).

[0132] In this invention, in the initial stage of iterative calculation, in order to accelerate the convergence speed of the reconstructed image, a step size s can be given. (n) (l) Multiply by a constant greater than 1, but verification is required.

[0133]

[0134] Where ∈ represents any small number greater than 0. If energy spectrum data is available... If the decrease condition (21) is not satisfied, the step size s is reduced (n) (l).

[0135] In summary, the steps of the multi-spectrum CT iterative reconstruction method with adaptive direction and step size and scatter correction provided by the application can be summarized as follows:

[0136]

[0137] Experimental results

[0138] In order to illustrate the technical effect of the method provided by the application, the experimental results are as follows.

[0139] The model of the experiment is a Shepp-Logan model, which is composed of a water base and a bone base. Figure 2 and Figure 3 The left ideal image of the water base is shown in the following figure. Figure 2 is the water base image result (gray window [0, 1.55]), Figure 3 is the bone base image result (gray window [0, 1.05]), Figure 4 is the pseudo single energy (50keV) image result (gray window [0.0113, 0.0352]).

[0140] The MCGPU program is used to simulate the generation of high and low energy projection images with scatter and noise, and then the substrate decomposition experiment is performed. The simulation parameters are as follows: SOD=500mm, SDD=1000mm, the detector matrix is 512x512, the pixel size is 0.8mm, the number of projection angles is 360, and 2^38 photons are used to simulate the projection image at each angle using two energy spectrum curves of 80kv and 125kv respectively. Among them, SOD represents the distance from the X-ray source S to the measured object, and SDD represents the distance from the X-ray source S to the measured object. The high and low energy data are incompatible, and the interval is 0.5 degrees.

[0141] The projection data is directly reconstructed using the substrate decomposition reconstruction algorithm, one is direct reconstruction, and the other is to add scatter removal in the iterative process, that is, to use the iterative algorithm ADS-SC with adaptive direction and step size and scatter correction proposed by the application. Figure 2 and Figure 3 The second column and the third column of and respectively represent the reconstruction results of the traditional substrate decomposition reconstruction method and the reconstruction results of the ADS-SC method of the application. From Figure 2 , Figure 3 It can be seen from the results of and that the image reconstructed by ADS-SC is closer to the ideal image.

[0142] The algorithm of substrate iteration is performed for 10 decomposition iterations, and the time is 21s. The method proposed in the application adds a scatter removal module on this basis, and the time for one round of scatter removal is 120s. The total time for two rounds of iteration is 282s. From the decomposition water-based results and the synthesized pseudo-single-energy image Figure 4 ), the reconstruction results after removing the scatter are closer to the original image, which shows that the multi-spectral CT iterative reconstruction method ADS-SC with scatter correction of adaptive direction and step length proposed in the application is effective.

[0143] Finally, it should be pointed out that: the above examples are only used to illustrate the technical solutions of the application, but not to limit them. It should be understood by those skilled in the art that the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the application.

Claims

1. A multi-spectral CT iterative reconstruction method with scatter correction and adaptive direction and step size, characterized in that, Comprising: Step 1, reconstructing a model of multi-energy CT projection with scatter correction by the following formula: （1） in: In the first Along the path under each energy spectrum Multicolor projection data of X-rays with scattering correction. ,express The base material along the ray path The distribution projection value, , Number of base materials Indicates the first Density distribution of the seed material in the tested object This represents the transpose of the vector; E represents the energy of the photon. and These are the lowest and highest energies of a photon, respectively. It is the first The normalized equivalent energy spectrum of an energy spectrum; , For the first The linear decay coefficient of the seed material; For the first energy spectrum The set of ray paths, The number of energy spectra; Step 2, solving the increment of the base material distribution projection value by using an iterative algorithm angle by angle, and updating the current base material image into a new base material image; comprising: Let the current basis material distribution projection value be , denote the value of the estimated basis material distribution projection after iterations, denote the number of iterations, initialized at ; all rays at the current angle , the following are calculated and the scattering intensity and the following are calculated and the following are calculated ; where represents the polychromatic projection data of the X-rays with scatter correction along the path at the th energy spectrum after iterations, represents the increment of the base material distribution projection value; Let , update the new base material distribution projection value , and back projection to get the first iteration after the base material image update.

Images