Phase contrast imaging simulation method, computer device and computer readable storage medium

Abstract

The phase contrast imaging simulation method comprises the steps that a light source, an imaged object, an object plane and an image plane are sequentially arranged in a simulation environment in the photon emission direction, and the object plane is divided into a plurality of object plane sub-areas; the interaction between the photons and the imaged object in the process of penetrating through the imaged object is simulated; for each object plane sub-region in the plurality of object plane sub-regions, determining the energy and phase distribution of a wave surface corresponding to the object plane sub-region based on the energy and phase of all photons penetrating through the imaged object to reach the object plane sub-region; based on the energy and phase distribution of each wave surface corresponding to each object plane sub-region, determining the energy distribution of all wave surfaces reaching the image plane; and generating a simulated image of the imaged object based on the energy distribution from the wave surfaces corresponding to the plurality of object plane sub-regions to the image plane. According to the invention, while phase contrast imaging is accurately simulated, time and computing resources required by simulation operation can be reduced.

Description

TECHNICAL FIELD

[0001] The present disclosure relates to the field of computer simulation, and more particularly, to a phase-contrast imaging simulation method, a computer device and a computer-readable storage medium. BACKGROUND

[0002] Phase-contrast imaging (PCI) is an imaging technique that extracts phase information generated by light passing through an object and converts the phase information into signal intensity of an image. In order to meet the needs of research and development, as well as cost and risk control, the phase-contrast imaging operation can usually be simulated in advance, and therefore, the phase-contrast imaging simulation technology is introduced into the fields of biomedical science and material science. However, in the actual application process, how to simulate high-resolution high-energy phase-contrast imaging under the condition of limited computing resources is still a problem to be solved in the industry. SUMMARY

[0003] In order to at least partially solve the above and other possible problems, embodiments of the present disclosure provide a phase-contrast imaging simulation scheme.

[0004] According to a first aspect of the present disclosure, a phase-contrast imaging simulation method is provided, which can include: sequentially arranging, in a simulation environment, a light source emitting photons, an imaged object, an object plane and an image plane along a photon emission direction, wherein the object plane is divided into a plurality of object plane sub-regions; simulating interaction of the photons with the imaged object in the process of passing through the imaged object; for each object plane sub-region of the plurality of object plane sub-regions, determining energy and phase distribution of a wave surface corresponding to the object plane sub-region based on energy and phase of all photons passing through the imaged object to reach the object plane sub-region; determining energy distribution of the wave surface corresponding to all object plane sub-regions in the object plane reaching the image plane based on the energy and phase distribution of each wave surface corresponding to each object plane sub-region; and generating simulated imaging of the imaged object based on the energy distribution of the wave surface corresponding to the plurality of object plane sub-regions reaching the image plane.

[0005] The present disclosure divides the object plane into a plurality of object plane sub-regions and regards each object plane sub-region as a corresponding wave surface to calculate the energy distribution of each wave surface reaching the image plane, thereby realizing the simulation operation of the phase-contrast imaging of the imaged object, and thus reducing the time and computing resources required for the simulation operation while accurately simulating the phase-contrast imaging.

[0006] In embodiments of the present disclosure, determining the energy and phase distribution of the wavefront corresponding to the object plane sub-region can include determining the energy of the wavefront based on the sum of the energies of all the photons reaching the object plane sub-region, and determining the phase of each position in the wavefront based on the average phase of all the photons reaching the object plane sub-region to obtain the phase distribution.

[0007] In embodiments of the present disclosure, determining the phase of each position in the wavefront can include, for each of the plurality of object plane sub-regions, determining a phase gradient of the wavefront corresponding to the object plane sub-region based on at least two object plane sub-regions close to the object plane sub-region and / or based on the object plane sub-region and at least one object plane sub-region close to the object plane sub-region, and determining the phase of each position in the wavefront based on the average phase of all the photons reaching the object plane sub-region and the phase gradient.

[0008] In embodiments of the present disclosure, determining the phase gradient of the wavefront corresponding to the object plane sub-region can include determining the phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between the average phase of each of two object plane sub-regions close to the object plane sub-region and a position difference between the geometric centers of each of the two object plane sub-regions.

[0009] In embodiments of the present disclosure, determining the phase gradient of the wavefront corresponding to the object plane sub-region can include determining the phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between the average phase of another object plane sub-region immediately adjacent to the object plane sub-region and the average phase of yet another object plane sub-region immediately adjacent to the another object plane sub-region and a position difference between the geometric centers of each of the another object plane sub-region and the yet another object plane sub-region.

[0010] In embodiments of the present disclosure, determining the phase gradient of the wavefront corresponding to the object plane sub-region can include determining the phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between the average phase of the object plane sub-region and the average phase of another object plane sub-region close to the object plane sub-region and a position difference between the geometric centers of each of the object plane sub-region and the another object plane sub-region.

[0011] In embodiments of the present disclosure, determining the phase of each position in the wave front based on the average phase of all the photons reaching the sub-region of the object plane and the phase gradient can include: taking the average phase as the phase of the geometric center position of the wave front; and for other positions in the wave front other than the geometric center position, determining the phase of each other position based on the position difference between each other position and the geometric center position, the phase gradient, and the phase of the geometric center position.

[0012] In embodiments of the present disclosure, determining the energy distribution can include: dividing the image plane into a plurality of image plane pixel regions; for each image plane pixel region in the plurality of image plane pixel regions, determining the wave field value of each wave front on the plurality of image plane pixel regions when the wave front reaches the image plane based on the energy and phase distribution of the wave front corresponding to each sub-region of the object plane; superimposing the wave field value of each wave front reaching each image plane pixel region to obtain the total wave field value for the image plane pixel region; and determining the total energy of the wave front corresponding to all sub-regions of the object plane reaching each image plane pixel region based on the total wave field value to determine the energy distribution, and the generating the simulated imaging of the imaged object can include: determining the gray value of each image plane pixel region based on the total energy corresponding to each image plane pixel region to generate the simulated imaging of the imaged object.

[0013] In embodiments of the present disclosure, the wave field value of each wave front on each of the plurality of image plane pixel regions when the wave front reaches the image plane can be determined by the following formula:

[0014]

[0015] wherein U(x, y, R2)|m, n represents the wave field value of the wave front corresponding to the mth row and nth column sub-region of the object plane on the image plane pixel region with position coordinates (x, y) at the geometric center position when the wave front reaches the image plane at a distance R2 from the object plane, R1 is the distance from the object plane to the light source, R2 is the distance from the object plane to the image plane, λ is the wavelength of the photon, k = 2π / λ, (m, n) is the number of the sub-region of the object plane, ξ m is the abscissa of the geometric center position of the mth row and nth column sub-region of the object plane, η n is the ordinate of the geometric center position of the mth row and nth column sub-region of the object plane, A(ξ m , η n) is the amplitude of the wavefront, the energy of the wavefront is equal to the square of the amplitude of the wavefront, x0, y0 are the horizontal and vertical coordinates of the initial position of the photon generated by the light source, u, v are the horizontal and vertical coordinates of any position of the wavefront, Φ m，n (u, v) is the phase of the (u, v) coordinate position of the wavefront, dξ is the distance between the geometric center positions of two adjacent sub-regions of the object plane in the horizontal coordinate, dη is the distance between the geometric center positions of two adjacent sub-regions of the object plane in the vertical coordinate, e represents a natural constant, and i represents an imaginary unit.

[0016] In an embodiment of the present disclosure, simulating the interaction of the photon with the imaged object during the process of passing through the imaged object includes simulating at least one of the following effects generated by the photon during the process of passing through the imaged object: coherent scattering, photoelectric effect, Compton effect, electron pair effect, refraction, and phase shift.

[0017] In an embodiment of the present disclosure, the object plane is arranged at a distance of 0 to 2 meters from the side of the imaged object close to the object plane.

[0018] In an embodiment of the present disclosure, the shapes and sizes of the divided plurality of sub-regions of the object plane are the same, and the shape of the sub-region of the object plane is one of the following: square, rectangle, triangle, hexagon, and circle.

[0019] According to a second aspect of the present disclosure, a computer device is provided, comprising a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method of the first aspect.

[0020] According to a third aspect of the present disclosure, a computer readable storage medium is provided, having a computer program stored thereon, wherein the computer program is executed by a processor to implement the steps of the method of the first aspect.

[0021] According to a fourth aspect of the present disclosure, a computer program product is provided, comprising a computer program, wherein the computer program is executed by a processor to implement the steps of the method of the first aspect.

[0022] The summary is provided to introduce a selection of concepts in a simplified form, which will be further described below in the detailed description. The summary is not intended to identify key or essential features of the present disclosure, nor is it intended to limit the scope of the present disclosure. BRIEF DESCRIPTION OF DRAWINGS

[0023] The above and other objects, features and advantages of the present disclosure will become more apparent from the following detailed description when taken in conjunction with the accompanying drawings in which like reference characters refer to like parts throughout the figures, and in which:

[0024] Figure 1 A flowchart of a phase contrast imaging simulation method according to an embodiment of the present disclosure is shown.

[0025] Figure 2 A schematic diagram of a simulation environment according to an embodiment of the present disclosure is shown.

[0026] Figure 3 A schematic diagram of an object plane according to an embodiment of the present disclosure is shown.

[0027] Figure 4 A schematic diagram of a portion of an object plane according to an embodiment of the present disclosure is shown.

[0028] Figure 5 A flowchart of a process for determining the phase of each location in a wavefront according to an embodiment of the present disclosure is shown.

[0029] Figure 6 A schematic diagram of a portion of an object plane used in determining the phase gradient of a wavefront according to an embodiment of the present disclosure is shown. DETAILED DESCRIPTION

[0030] Embodiments of the present disclosure will be described more fully hereinafter with reference to the accompanying drawings, in which embodiments of the present disclosure are shown. This disclosure may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art. One skilled in the art can readily devise many alternative technology solutions that fall within the scope of the present disclosure without departing from the spirit and scope of the present disclosure.

[0031] As used herein, the term "includes" and its variants are meant to be construed as encompassing, but not limited to, the following: "including, but not limited to"; "including and including, but not limited to"; "including, but not limited to and including"; "including other items in addition to those listed"; "including at least the items listed after this term"; "including, but limited to"; and "including at least the following items, including collectively those items listed after this term, other items in addition to those listed, and other items left open." Unless specifically stated otherwise, the term "or" as used herein is intended to mean "and / or," i.e., the term "or" as used herein is intended to be interpreted in the same way as the term "and / or" is interpreted. Unless specifically stated otherwise, the term "based on" as used herein is intended to mean "based, at least in part, on" and "based, at least in part, on." The terms "one example embodiment," "an example embodiment," "one embodiment," and "an embodiment" as used herein indicate the existence of at least one example embodiment and do not impose a limitation on the present disclosure, except insofar as a claim might limit a particular embodiment. Other embodiments of the present disclosure will be apparent to those skilled in the art from consideration of the specification and practice of the present disclosure disclosed herein. Therefore, it is intended that the specification and examples be considered as exemplary only with the true scope and spirit of the present disclosure being indicated by the following claims.

[0032] Further, as used herein, the term "adjacent" means that one object is in direct contact with another object, and the term "near" means that one object is disposed around another object and not in direct contact. The specific meanings of these two terms will be further explained below in conjunction with the drawings.

[0033] To simulate high-resolution high-energy phase contrast imaging under the condition of limited computing resources, the present disclosure provides a phase contrast imaging simulation method. Figure 1 A flowchart of the phase contrast imaging simulation method 100 according to an embodiment of the present disclosure is shown. It should be understood that, Figure 1 Any step in the flowchart can be executed in the form of a computer program on a computing device and stored on a storage device. The computing device can be any type of fixed computing device, mobile computing device, or portable computing device, including but not limited to a desktop computer, a laptop computer, a notebook computer, a netbook computer, a tablet computer, a mobile phone, etc. In addition, the phase contrast imaging mentioned in the present disclosure can be based on the propagation of phase contrast imaging, or other types of phase contrast imaging.

[0034] At S110, the computing device can sequentially arrange the light source, the imaged object, the object plane, and the image plane along the direction of photon emission in the simulation environment. In the simulation environment, the computing device can simulate the emission of photons by the light source, can simulate the interaction of the photons with the imaged object, and can simulate the energy and phase distribution of the photons transported to the object plane and the image plane.

[0035] To clearly describe the embodiments of the present disclosure, Figure 2 A schematic diagram of a simulation environment 200 according to an embodiment of the present disclosure is shown. As Figure 2 shown, the simulation environment 200 includes a light source 210 emitting photons, an imaged object 220, an object plane 230, and an image plane 240. These virtual components are sequentially arranged along the optical axis A (i.e., the direction of photon emission) of the light source 210. In Figure 2 , R1 represents the distance from the object plane 230 to the light source 210, R2 represents the distance from the object plane 230 to the image plane 240, the object plane 230 is perpendicular to the optical axis A, and the image plane 240 is also perpendicular to the optical axis A.

[0036] In some embodiments, the object plane 230 is at a distance R3 from the side of the imaged object 220 close to the object plane 230, and R3 is in the range of 0 to 2 meters, for example, R3 can be equal to 0 centimeters, i.e., the object plane 230 is immediately adjacent to the side of the imaged object 220 away from the light source 210. Alternatively or additionally, R3 is in the range of 0 to 1 centimeter, or in the range of 0 to 10 centimeters, or in the range of 0 to 50 centimeters, or in the range of 0 to 1 meter, or in the range of 0 to 3 meters.

[0037] To avoid the massive computation caused by calculating the transport process of each simulated photon one by one, the photons on each region can be analyzed locally on the object plane 230. In some embodiments, the object plane 230 can be pre-divided into multiple object plane sub-regions. Figure 3A schematic diagram of the object plane 230 according to an embodiment of the present disclosure is shown. It should be understood that... Figure 3 The object plane 230 in the middle can be considered as being from... Figure 2 The light source 210 or the plane observed by the imaged object 220 towards the object plane 230.

[0038] like Figure 3 As shown, object plane 230 comprises 16 object plane sub-regions arranged in 4 rows and 4 columns, such as object plane sub-region 2301 located in the upper right corner of object plane 230. It should be understood that object plane 230 can contain any number of object plane sub-regions, and the object plane sub-regions in object plane 230 can be arranged in any manner, and are not limited to... Figure 3 The objects shown are arranged in squares or rectangles. For the sake of convenience in describing embodiments of this disclosure, the object plane 230 can be considered to contain N. o ×N o N object plane sub-regions o It is a natural number. In some embodiments, the multiple object plane sub-regions divided in the object plane 230 have the same shape and size. It should be understood that the shapes and sizes of the multiple object plane sub-regions can also be set to be different. Furthermore, the shape of each object plane sub-region can be a square, rectangle, triangle, hexagon, circle, or other shapes. Object plane sub-regions with a square or rectangular shape can also be referred to as "grids" in this disclosure. Figure 3 In this simulation, the object plane 230 is divided into 4×4 grids. In this way, the energy and phase of photons that pass through the imaged object 220 and reach each grid can be statistically analyzed separately.

[0039] After setting up the light source 210, the object to be imaged 220, the object plane 230, and the image plane 240, the computing device can simulate the emission of photons by the light source 210. As an example, the computing device can randomly generate photons in a Monte Carlo simulation program based on the characteristics of the photons emitted by the light source 210 (e.g., energy spectrum, emission direction, etc.) input by the user. Furthermore, the object to be imaged 220 can be biological tissue, lesions, or other objects to be observed in medical research; therefore, the Monte Carlo simulation program also requires the user to input optical information about the object 220, such as its refractive index. It should be understood that, in addition to the Monte Carlo simulation program, other simulation programs that have been developed or will be developed can also simulate the operation of the light source 210 emitting photons.

[0040] Further, at S120, the computing device can simulate the interaction between the photons and the imaged object 220 during the process of the photons passing through the imaged object 220. It should be understood that the interaction between the photons and the imaged object 220 is complicated and random, which makes the Monte Carlo simulation program an effective statistical tool. For example, the Monte Carlo simulation program can simulate the interaction between the photons and the imaged object 220 and the motion trajectory of the photons after passing through the imaged object 220 according to the characteristics of the photons and the optical information of the imaged object 220.

[0041] As an example, after the Monte Carlo simulation program randomly generates a large number of photons, each photon is assigned an initial state, including the position, direction of travel, energy, phase, and the like of the photon. Then, the Monte Carlo simulation program calculates and updates the next state of the photon according to the optical information of the imaged object 220 and the initial state of the photon.

[0042] During the transport of the photons, if the photons pass through the imaged object 220, refraction, scattering, diffraction, and absorption may occur (multiple occurrences of the above phenomena may occur, for example, the photon is first refracted and then scattered, or the photon may eventually disappear due to absorption). At this time, the Monte Carlo simulation program calls the nuclear physics database to query the probabilities of coherent scattering effects, photoelectric effects, Compton effects, electron pair effects, refraction, and phase shift effects. Thus, the Monte Carlo simulation program can simulate the coherent scattering effects, photoelectric effects, Compton effects, electron pair effects, refraction, and phase shift effects generated by the photons during the process of passing through the imaged object 220 according to the data queried in the nuclear physics database, further simulate the transport process of the photons, and finally record the photons reaching the object plane 230 after the imaged object 220. It should be understood that in addition to the Monte Carlo simulation program, other simulation programs that have been developed or will be developed can also simulate the interaction between the light source 210 and the imaged object 220.

[0043] At S130, for each of the plurality of object plane sub-regions, the computing device can determine the energy and phase distribution of the wave surface corresponding to the object plane sub-region based on the energy and phase of all photons passing through the imaged object 220 and reaching the object plane sub-region. As an example, the computing device counts the energy and phase of all photons received by each object plane sub-region, and then for an object plane sub-region, the computing device regards all photons received by the object plane sub-region as a whole as a wave surface, so that the subsequent transport process of each wave surface can be further analyzed. Therefore, compared with the traditional simulation method of analyzing the transport process of each photon one by one, the simulation method of the present disclosure can save a large amount of computing resources and reduce the operation time by analyzing the transport process of the wave surface.

[0044] Figure 4A schematic diagram of a partial object plane sub-region is shown according to an embodiment of the present disclosure. As Figure 4 shown, the Monte Carlo simulation program records the photons that pass through the imaged object 220 to arrive on one object plane sub-region 410 in the object plane. In Figure 4 the object plane sub-region 410 corresponds to one wavefront, in other words, the set of all photons in the object plane sub-region 410 can be considered as one wavefront, and these photons and the data recorded by the Monte Carlo program about these photons such as energy, phase, etc. will be used to further calculate the energy and phase distribution of the wavefront. Then, each wavefront corresponding to each object plane sub-region respectively is considered as a whole, so that the transport of each wavefront towards the image plane 240 can be analyzed and simulated. It should be understood that, Figure 4 the black dots in the above are only used to exemplarily show the photons that arrive at the object plane sub-region 410, and are not intended to show the number, energy, phase, or even transport direction of the photons.

[0045] In some embodiments, the computing device can determine the sum of the energy of all the photons that arrive at the object plane sub-region 410 as the energy of the wavefront corresponding to the object plane sub-region 410. In other words, all the photons in the object plane sub-region 410 will be propagated as a whole to the image plane 240 as one wavefront, and the energy of any position in the wavefront is equivalent to the energy of the wavefront, i.e. the sum of the energy of all the photons that arrive at the object plane sub-region 410. In addition, the energy of the wavefront is equal to the square of the amplitude of the wavefront. Therefore, after the energy of each wavefront is determined, the amplitude of each wavefront can be determined for the subsequent analytical calculation of the Fresnel diffraction integral.

[0046] In addition to determining the energy of the wavefront, the computing device can also determine the phase of each position in the wavefront based on the average phase of all the photons that arrive at the object plane sub-region 410 to obtain a phase distribution. Figure 5 A flowchart of a process 500 for determining the phase of each position in the wavefront according to an embodiment of the present disclosure is shown.

[0047] In S510, for each object plane sub-region of the plurality of object plane sub-regions, the computing device can determine the phase gradient of the wavefront corresponding to the object plane sub-region based on at least two object plane sub-regions close to the object plane sub-region and / or based on the object plane sub-region and at least one object plane sub-region close to the object plane sub-region. In order to facilitate understanding of the embodiments of the present disclosure, Figure 6 A schematic diagram of the object plane sub-regions used when determining the phase gradient of the wavefront according to an embodiment of the present disclosure is shown.

[0048] As Figure 6In some embodiments, the phase gradient of the wavefront corresponding to a certain object plane sub-region 510 is determined based on at least two object plane sub-regions close to the object plane sub-region 510. In this regard, the computing device can determine two object plane sub-regions close to the object plane sub-region 510 and determine the average phase of each of the two object plane sub-regions and the geometric center of each of the two object plane sub-regions. Then, the computing device can determine the phase gradient of the wavefront corresponding to the object plane sub-region 510 based on the difference between the two average phases and the difference in position between the two geometric centers. It can be appreciated that the closer the two selected object plane sub-regions are to the object plane sub-region 510, the closer the phase gradient calculated using the two selected object plane sub-regions can be to the phase gradient between different positions in the object plane sub-region 510, and thus the more accurate the phase values of different positions in the object plane sub-region 510 calculated using the phase gradient can be.

[0049] In some embodiments, the computing device can determine two object plane sub-regions immediately adjacent to the object plane sub-region 510, such as the object plane sub-regions 520, 530, or the object plane sub-regions 550, 560, or 530, 560, etc., to calculate the phase gradient for the object plane sub-region 510. In some embodiments, the computing device can also determine one object plane sub-region immediately adjacent to the object plane sub-region 510 and another object plane sub-region immediately adjacent to the one object plane sub-region, such as the object plane sub-regions 530, 540, to calculate the phase gradient for the object plane sub-region 510. In some embodiments, the computing device can also determine two object plane sub-regions located around the object plane sub-region 510, such as the object plane sub-regions 540, 570, etc., to calculate the phase gradient for the object plane sub-region 510.

[0050] In some embodiments, the phase gradient of the wavefront corresponding to a certain object plane sub-region is determined based on the object plane sub-region and at least one object plane sub-region close to the object plane sub-region. In this regard, the computing device can determine the object plane sub-region 510 and another object plane sub-region close to the object plane sub-region 510 and determine the average phase of each of the two object plane sub-regions and the geometric center of each of the two object plane sub-regions. Then, the computing device can determine the phase gradient of the wavefront corresponding to the object plane sub-region 510 based on the difference between the two average phases and the difference in position between the two geometric centers.

[0051] In some embodiments, the computing device can determine the object plane sub-region 510 and one object plane sub-region immediately adjacent to the object plane sub-region 510, e.g., the object plane sub-region 530, or the object plane sub-region 560, etc., to calculate the phase gradient for the object plane sub-region 510. In some embodiments, the computing device can determine the object plane sub-region 510 and one object plane sub-region located around the object plane sub-region 510, e.g., the object plane sub-region 570, or the object plane sub-region 540, etc., to calculate the phase gradient for the object plane sub-region 510.

[0052] In some embodiments, the lateral phase gradient and the longitudinal phase gradient in the coordinate system of the object plane can be calculated for a certain object plane sub-region 510, respectively, to calculate the phase of the wave front at an arbitrary position in the object plane sub-region 510 relative to the geometric center position of the object plane sub-region 510 based on the lateral phase gradient and the longitudinal phase gradient.

[0053] As an example, the phase gradient of the wave front in the lateral direction can be determined by using the following formula

[0054]

[0055] And the phase gradient of the wave front in the longitudinal direction can be determined by using the following formula.

[0056]

[0057] wherein, is the phase gradient of the wave front in the lateral direction, is the phase gradient of the wave front in the longitudinal direction, m, n are the numbers of the object plane sub-regions, which can represent the object plane sub-region in the mth row and the nth column, e.g., 0≤m o and 0≤n o . dξ is the distance between the geometric center positions of two object plane sub-regions immediately adjacent in the lateral coordinate axis direction of the object plane coordinate system, and dη is the distance between the geometric center positions of two object plane sub-regions immediately adjacent in the longitudinal coordinate axis direction of the object plane coordinate system. Therefore, the phase gradient of the wave front in the lateral coordinate axis can be determined by the difference between the average phases of the object plane sub-regions 520, 530 immediately adjacent to the left and right sides of the object plane sub-region 510, respectively, and the position difference between the geometric centers of the two object plane sub-regions, and the phase gradient of the wave front in the longitudinal coordinate can be determined by the difference between the average phases of the object plane sub-regions 550, 560 immediately adjacent to the upper and lower sides of the object plane sub-region 510, respectively, and the position difference between the geometric centers of the two object plane sub-regions. In this way, the phase of each position of the wave front can be further obtained by linear fitting, etc.

[0058] Afterwards, in S520, the computing device can determine the phase of each position in the wavefront corresponding to the object plane sub-region 510 based on the average phase of all the photons reaching the object plane sub-region 510 and the phase gradient obtained above. In some embodiments, the computing device can take the average phase of the object plane sub-region 510 as the phase of the geometric center position of the wavefront corresponding to this object plane sub-region 510, and for other positions in the wavefront other than the geometric center position, the computing device can determine the phase of each other position based on the position difference between each other position and the geometric center position, the phase gradient obtained above, and the phase of the geometric center position.

[0059] As an example, the phase of other positions on the wavefront corresponding to the object plane sub-region 510 other than the geometric center position can be determined by linear fitting using the following formula.

[0060]

[0061] where m, n are the indices of the wavefront corresponding to the object plane sub-region as described above, ξ is the horizontal coordinate of the other position on the wavefront other than the geometric center position, η is the vertical coordinate of the other position, Φ m,n (ξ, η) is the phase at the position with coordinates (ξ, η) on the wavefront corresponding to the object plane sub-region in the mth row and nth column of the object plane 230. Φ m,n is the phase of the center position of the wavefront, ξ m is the horizontal coordinate of the geometric center position of the wavefront, η n is the vertical coordinate of the geometric center position of the wavefront.

[0062] The above embodiments determine the phase distribution of the wavefront corresponding to the object plane sub-region in a fitting manner by determining the average phase of each wavefront as the phase of the center position of the wavefront and determining the phase gradient using adjacent wavefronts. This process can directly calculate the phase of any position on the wavefront using formula 3, which reduces the computational complexity relative to the traditional way of directly performing integration operations. In addition, the scheme of the present disclosure can reduce the number or density of object plane sub-regions, saving computing resources. Furthermore, the above embodiments directly determine the sum of the energies of all the photons reaching an object plane sub-region as the energy of the wavefront corresponding to the object plane sub-region, thereby simplifying the energy calculation.

[0063] Afterwards, returning to Figure 1 At S140, the computing device can determine the energy distribution of the wavefront corresponding to all the object plane sub-regions in the object plane 230 reaching the image plane 240 based on the energy and phase distribution of each wavefront corresponding to each object plane sub-region.

[0064] Specifically, the computing device can divide the image plane 240 into multiple image plane pixel regions. It should be understood that dividing the image plane 240 can be done in a similar manner to dividing the object plane 230. It is important to note that while the size of the object plane sub-regions can be arbitrarily changed within a certain range according to user needs when dividing the object plane 230, the size of the image plane pixel regions typically depends on the required image pixel size. Therefore, the size and number of the divided image plane pixel regions can be determined based on the area of ​​the image plane and the required imaging resolution.

[0065] In order to determine the energy distribution of all wavefronts reaching the image plane 240, the computing device needs to determine the wavefield value of each wavefront in each of the multiple image plane pixel regions that are divided into multiple image plane pixel regions when it reaches the image plane.

[0066] As an example, we can substitute Equation 3 into the Fresnel diffraction integral formula to obtain the following formula:

[0067]

[0068] Where, U(x, y, R2)| m，n This can represent the wavefield value on the image plane pixel region with position coordinates (x, y) at the geometric center when the wavefront corresponding to the m-th row and n-th column sub-region of object plane 230 reaches the image plane 240 at a distance R2 from that object plane. The wavefield value is a physical quantity describing the amplitude, phase, and other characteristics of the wavefront at a certain position in space. R1 is the distance from object plane 230 to light source 210, R2 is the distance from object plane 230 to image plane 240, λ is the wavelength of the photon, k is the wave number, k = 2π / λ, (m, n) is the number of the object plane sub-region, ξ m Let η be the x-coordinate of the geometric center of the object plane sub-region in the m-th row and n-th column. n Let A(ξ) be the ordinate of the geometric center position of the object plane sub-region in the m-th row and n-th column. m η n Let y0 be the amplitude of the wavefront, where the energy of the wavefront is equal to the square of the amplitude, x0 and y0 are the initial x and y coordinates of the photon generated by the light source 210, u and v are the integration variables, representing the x and y coordinates of any position on the wavefront, and Φ is the amplitude of the wavefront. m，n (u, v) represents the phase of the wavefront at the (u, v) coordinate position, dξ represents the distance between the geometric centers of two adjacent object plane sub-regions on the horizontal axis, dη represents the distance between the geometric centers of two adjacent object plane sub-regions on the vertical axis, e represents the natural constant, and i represents the imaginary unit.

[0069] It should be understood that (x, y) is not limited to representing the position coordinates of the geometric center position of the image plane pixel region, but can also be the position coordinates of the edge position (for example, the lower left corner) of the image plane pixel region.

[0070] In addition, since formula 4 is a Fresnel diffraction integral that directly uses formula 3 to calculate the wave surface phase, formula 4 can be further analytically calculated, so that its analytical expression is as follows:

[0071]

[0072] wherein, and, Erfi() is the imaginary error function in the complex domain, and z and z0 are specific position coordinate variables, and in formula 5, x and x0 are respectively brought into z and z0 for operation. It should be understood that the traditional Fresnel diffraction integral considering the phase of any position on the wave surface needs to integrate the phase therein, so the traditional Fresnel diffraction integral needs considerable amount of operation and has high calculation complexity. The present disclosure directly determines the phase of any position on each wave surface by using formula 3, and then substitutes formula 3 into formula 4, so as to obtain the analytical expression formula 5 of the Fresnel diffraction integral. Therefore, the present disclosure can solve the Fresnel diffraction integral based on the algebraic operation of formula 5, which avoids the complexity of calculation and reduces the amount of calculation.

[0073] After the wave field value is determined, the computing device can further superimpose the wave field value of each wave surface reaching each image plane pixel region to obtain the total wave field value for the image plane pixel region by the following formula.

[0074]

[0075] wherein S(x, y, R2) can represent the superposition of the wave field values of the wave surfaces corresponding to all N o ×N o sub-regions of the object plane 230 on the image plane pixel region with the position coordinates (x, y) of the geometric center position when reaching the image plane 240 at a distance R2 from the object plane. m，n U(x, y, R2) can represent the wave field value of the wave surface corresponding to the mth row and nth column sub-region of the object plane 230 on the image plane pixel region with the position coordinates (x, y) of the geometric center position when reaching the image plane 240 at a distance R2 from the object plane.

[0076] Afterwards, the computing device can determine the total energy of the wavefronts corresponding to all the sub-regions of the object plane 230 reaching each of the pixel regions of the image plane 240 based on the total wavefield values solved above to determine the energy distribution of the wavefronts corresponding to all the sub-regions of the object plane 230 reaching the image plane 240. In some embodiments, the total energy of the wavefronts corresponding to all the sub-regions of the object plane 230 reaching each of the pixel regions of the image plane 240 can be determined based on the following equation.

[0077] I(x, y, R2) = [S(x, y, R2)] 2 Equation 7

[0078] where I(x, y, R2) can represent the total energy of the wavefronts corresponding to all the sub-regions of the object plane 230 reaching the image plane 240 at a distance R2 from the object plane at the geometric center position of the pixel region of the image plane with the position coordinates (x, y).

[0079] It can be understood that other operation results associated with the total wavefield values solved above can also be determined as the total energy described above. As an example, after the total wavefield values corresponding to each of the pixel regions of the image plane are determined, the total wavefield values can be directly used as the total energy for subsequent calculations. As another example, a value that is a predetermined multiple of the total wavefield values can also be determined as the total energy described above.

[0080] At S150, the computing device can generate a simulated imaging of the imaged object 220 based on the energy distribution of the wavefronts corresponding to all the sub-regions of the object plane 230 reaching the image plane 240. In some embodiments, to generate the simulated imaging of the imaged object 220, the computing device can determine the gray scale values of each of the pixel regions of the image plane based on the total energy corresponding to each of the pixel regions of the image plane to generate the simulated imaging of the imaged object 220. Alternatively or additionally, the signal strength or brightness values of each of the pixel regions of the image plane can also be determined based on the total energy determined above. Since each of the pixels of the image plane 240 has a determined corresponding gray scale or brightness, the simulated imaging can be generated.

[0081] In addition, in some embodiments, the user can pre-set the number of photons to be simulated, and when the number of photons participating in the simulation is the same as the pre-set number, the computing device can generate a plot based on the simulated imaging and show it to the user.

[0082] In some embodiments, to make the simulated imaging closer to the actual result, various noise and image blur factors can be further added to blur the ideal energy distribution solved above to obtain a simulated result that is closer to the real result.

[0083] The phase contrast imaging simulation method proposed in the present disclosure can solve the problems of high computational complexity and high simulation overhead faced in simulating high-resolution high-energy phase contrast imaging. The core of the method proposed in the present disclosure is to divide the object plane into multiple sub-regions and regard all photons on each sub-region as a corresponding wave front, thereby converting the simulation of the transport of multiple photons from the object plane to the image plane in the particle dimension into simulation in the wave front dimension. In addition, during the simulation, the present disclosure determines the total energy of all photons of each object plane sub-region as the energy of the corresponding wave front, and determines the phase distribution of each position within each wave front by linear fitting based on the average phase determined by the phase of all photons of each object plane sub-region, so that the analytical wave front Fresnel diffraction integral can be introduced, so as to include the parameters related to the photon energy in the analytical integral, so that the numerical calculation is decoupled from the energy. This is equivalent to replacing the numerical integral operation with an algebraic operation, thereby significantly improving the computational efficiency and saving the computing resources. In addition, compared with the traditional Monte Carlo simulation method, this method considers the interference effect through the calculation of the Fresnel diffraction integral, so that the simulation result is more realistic. In summary, the phase contrast imaging simulation method of the present disclosure improves the realism of the simulation result while reducing the computational complexity and cost.

[0084] It should be understood that, although several means or sub-means are mentioned in the foregoing detailed description of the apparatus, such division is merely exemplary and not mandatory. Indeed, according to embodiments of the present disclosure, features and functions of two or more means described above can be embodied in one means. Conversely, features and functions of one means described above can be further divided into means embodied by several means.

[0085] Many modifications and other embodiments of the present disclosure set forth herein will be apparent to those of ordinary skill in the art from the teachings herein. Thus, it is intended that the embodiments of the present disclosure not be limited to the particular embodiments disclosed, and that modifications and other embodiments are intended to fall within the scope of the present disclosure. Further, while the above description has been made based on certain example combinations of components and / or functions, it is to be appreciated that different combinations of components and / or functions can be provided by alternative embodiments without departing from the scope of the present disclosure. In this regard, for example, other combinations of components and / or functions differing from those explicitly described above are also contemplated as within the scope of the present disclosure. While specific terminology has been employed herein, such terminology is used in the generic and descriptive sense only, and not for purposes of limitation.

Claims

1. A method of phase-contrast imaging simulation, comprising: sequentially arranging, in a simulation environment, a light source emitting photons, an imaged object, an object plane and an image plane along a direction of photon emission, wherein the object plane is divided into a plurality of object plane sub-regions; simulating interactions of the photons with the imaged object in a process of passing through the imaged object; for each of the plurality of object plane sub-regions, determining an energy and phase distribution of a wavefront corresponding to the object plane sub-region based on energies and phases of all photons passing through the imaged object to reach the object plane sub-region; determining an energy distribution of the wavefronts corresponding to all of the object plane sub-regions of the object plane reaching the image plane based on the energy and phase distribution of each wavefront corresponding to each of the object plane sub-regions; and generating a simulated imaging of the imaged object based on the energy distribution of the wavefronts corresponding to the plurality of object plane sub-regions reaching the image plane.

2. The method of phase-contrast imaging simulation of claim 1, wherein determining the energy and phase distribution of the wavefront corresponding to the object plane sub-region comprises: determining the energy of the wavefront based on a sum of energies of all photons reaching the object plane sub-region; and determining a phase of each location in the wavefront based on an average phase of all photons reaching the object plane sub-region to obtain a phase distribution.

3. The method of phase-contrast imaging simulation of claim 2, wherein determining the phase of each location in the wavefront comprises: for each of the plurality of object plane sub-regions, determining a phase gradient of the wavefront corresponding to the object plane sub-region based on at least two object plane sub-regions close to the object plane sub-region and / or based on the object plane sub-region and at least one object plane sub-region close to the object plane sub-region; and determining the phase of each location in the wavefront based on the average phase of all photons reaching the object plane sub-region and the phase gradient.

4. The method of phase-contrast imaging simulation of claim 3, wherein determining the phase gradient of the wavefront corresponding to the object plane sub-region comprises: determining the phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between respective average phases of two object plane sub-regions close to the object plane sub-region and a position difference between respective geometric centers of the two object plane sub-regions.

5. The method of phase-contrast imaging simulation of claim 4, wherein determining the phase gradient of the wavefront corresponding to the object plane sub-region comprises: determining the phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between an average phase of another object plane sub-region immediately adjacent to the object plane sub-region and an average phase of yet another object plane sub-region immediately adjacent to the another object plane sub-region and a position difference between respective geometric centers of the another object plane sub-region and the yet another object plane sub-region.

6. The method of phase-contrast imaging simulation of claim 3, wherein determining the phase gradient of the wavefront corresponding to the object plane sub-region comprises: ​ ​ determining a phase gradient of the wavefront corresponding to the object plane sub-region based on a difference between an average phase of the object plane sub-region and an average phase of another object plane sub-region close to the object plane sub-region and a position difference between a geometric center of each of the object plane sub-region and the another object plane sub-region.

7. The method of claim 3, wherein determining the phase of each location in the wavefront based on the average phase and the phase gradient of all the photons reaching the object plane sub-region comprises: using the average phase as the phase of a geometric center location of the wavefront; for other locations in the wavefront other than the geometric center location, determining the phase of each other location based on a position difference between each other location and the geometric center location, the phase gradient, and the phase of the geometric center location.

8. The method of claim 1, wherein determining the energy distribution comprises: dividing the image plane into a plurality of image plane pixel regions; for each image plane pixel region of the plurality of image plane pixel regions, determining wavefield values of each wavefront reaching the plurality of image plane pixel regions based on the energy and phase distribution of the wavefront corresponding to each object plane sub-region of the plurality of object plane sub-regions; superimposing the wavefield values of each wavefront reaching each of the image plane pixel regions to obtain a total wavefield value for the image plane pixel region; and determining a total energy of the wavefront corresponding to all of the object plane sub-regions of the object plane reaching each of the image plane pixel regions based on the total wavefield values to determine the energy distribution, and, the generating the simulated imaging of the imaged object comprises: determining a gray value of each of the image plane pixel regions based on the total energy corresponding to each of the image plane pixel regions to generate the simulated imaging of the imaged object.

9. The method of claim 8, wherein the wavefield values of each wavefront reaching each of the plurality of image plane pixel regions are determined by the following equation:

10. The method of claim 1, wherein simulating the interaction of the photons with the imaged object during the process of passing through the imaged object comprises simulating at least one of the following effects produced by the photons during the process of passing through the imaged object: wherein U(x, y, R2) m,n where U(x, y, R2) represents the wave field value on the image plane pixel region with position coordinate (x, y) at the geometric center position when the wave front corresponding to the mth row and nth column object plane sub-region of the object plane reaches the image plane with a distance R2 from the object plane R2, R1 is the distance from the object plane to the light source, R2 is the distance from the object plane to the image plane, λ is the wavelength of the light particle, k = 2π / λ, (m, n) is the number of the object plane sub-region, ξ m is the horizontal coordinate of the geometric center position of the mth row and nth column object plane sub-region, η n is the vertical coordinate of the geometric center position of the mth row and nth column object plane sub-region, A(ξ m , η n ) is the amplitude of the wave front, the energy of the wave front is equal to the square of the amplitude of the wave front, x0, y0 are the horizontal and vertical coordinates of the initial position of the light particle generated by the light source, u, v are integral variables, φ m，n (u, v) is the phase of the (u, v) coordinate position of the wave front, dξ is the distance between the geometric center positions of two adjacent object plane sub-regions in the horizontal coordinate, dη is the distance between the geometric center positions of two adjacent object plane sub-regions in the vertical coordinate, e represents the natural constant, and i represents the imaginary unit. coherent scattering; photoelectric effect; Compton effect; pair production; refraction; and phase shift.

11. The method of claim 1, wherein the object plane is set to be 0 to 2 meters away from a side of the imaged object close to the object plane.

12. The method of claim 8, wherein the plurality of object plane sub-regions are of the same shape and size, and wherein the shape of the object plane sub-region is one of: square, rectangle, triangle, hexagon, circle. the processor executes the computer program to implement the steps of the method of any one of claims 1-12.

13. A computer apparatus comprising a memory, a processor, and a computer program stored on the memory, wherein, 13. A computer program product comprising a computer readable medium storing instructions that, when executed by a processor, cause the processor to perform the steps of the method of any one of claims 1-12.

14. A computer readable storage medium having stored thereon a computer program, wherein, The computer program, which when executed by the processor, implements the steps of the method of any one of claims 1-12.

15. A computer program product comprising a computer program, wherein, The computer program, which when executed by the processor, implements the steps of the method of any one of claims 1-12.

Images