CN113284205B - CT iterative reconstruction method and device - Google Patents

CT iterative reconstruction method and device Download PDF

Info

Publication number
CN113284205B
CN113284205B CN202110443876.3A CN202110443876A CN113284205B CN 113284205 B CN113284205 B CN 113284205B CN 202110443876 A CN202110443876 A CN 202110443876A CN 113284205 B CN113284205 B CN 113284205B
Authority
CN
China
Prior art keywords
integral
value
image
weight matrix
image sequence
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.)
Active
Application number
CN202110443876.3A
Other languages
Chinese (zh)
Other versions
CN113284205A (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.)
Neusoft Medical Systems Co Ltd
Original Assignee
Shenyang Advanced Medical Equipment Technology Incubation Center Co ltd
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 Shenyang Advanced Medical Equipment Technology Incubation Center Co ltd filed Critical Shenyang Advanced Medical Equipment Technology Incubation Center Co ltd
Priority to CN202110443876.3A priority Critical patent/CN113284205B/en
Publication of CN113284205A publication Critical patent/CN113284205A/en
Application granted granted Critical
Publication of CN113284205B publication Critical patent/CN113284205B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration by the use of more than one image, e.g. averaging, subtraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10016Video; Image sequence
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20212Image combination
    • G06T2207/20224Image subtraction

Abstract

The embodiment of the invention provides a CT iterative reconstruction method and device. In the embodiment of the invention, the three-dimensional weight matrix required by three-dimensional forward projection is regarded as the product of an XY plane weight matrix and a Z-direction weight vector, and the XY plane weight matrix is determined and stored in advance in an offline manner, and the method comprises the following steps of: reading a stored XY plane weight matrix, performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data, determining an output image sequence after the current iteration according to the projection data, taking the output image sequence after the current iteration as a final CT reconstruction image sequence if a preset iteration stop condition is met, otherwise, taking the output image sequence after the current iteration as an input image sequence of the next iteration, executing the next iteration, and replacing the generation operation of the matrix by using the operation of reading the XY plane weight matrix in the iteration process, thereby reducing iteration time consumption and improving the reconstruction speed.

Description

CT iterative reconstruction method and device
Technical Field
The invention relates to the technical field of medical image processing, in particular to a CT iterative reconstruction method and device.
Background
CT (Computed Tomography), computed tomography) imaging is a very widely used imaging technique. The CT analysis reconstruction algorithm and the CT iteration reconstruction algorithm are two common CT image reconstruction modes.
Compared with a CT analytic reconstruction algorithm, the CT iterative reconstruction algorithm has the advantages that a CT imaging system can be modeled by combining prior information, and the CT analytic reconstruction algorithm is suitable for the conditions of insufficient information quantity and high noise. The disadvantage of the CT iterative reconstruction algorithm is that it is time consuming to calculate and slow.
Disclosure of Invention
In order to overcome the problems in the related art, the invention provides a CT iterative reconstruction method and device, a CT device and a CT system, and the reconstruction speed is improved.
According to a first aspect of the embodiment of the present invention, a CT iterative reconstruction method is provided, wherein a three-dimensional weight matrix required by three-dimensional forward projection is regarded as a product of an XY plane weight matrix and a Z direction weight vector, and the XY plane weight matrix is determined offline in advance and stored;
the iterative reconstruction process comprises the following steps:
reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
determining an output image sequence after the iteration according to the projection data;
if the preset iteration stopping condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
According to a second aspect of an embodiment of the present invention, there is provided a CT iterative reconstruction apparatus regarding a three-dimensional weight matrix required for three-dimensional forward projection as a product of an XY plane weight matrix and a Z direction weight vector, the apparatus storing therein the XY plane weight matrix determined offline in advance, the apparatus comprising:
the projection module is used for reading the stored XY plane weight matrix in each iterative reconstruction process, and carrying out three-dimensional forward projection on the input image sequence by utilizing the XY plane weight matrix to obtain projection data;
the determining module is used for determining an output image sequence after the iteration according to the projection data;
the output module is used for taking the output image sequence after the iteration as a final CT reconstructed image sequence if a preset iteration stop condition is met; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
The technical scheme provided by the embodiment of the invention can have the following beneficial effects:
in the embodiment of the invention, the three-dimensional weight matrix required by three-dimensional forward projection is regarded as the product of an XY plane weight matrix and a Z-direction weight vector, and the XY plane weight matrix is determined and stored in advance in an offline manner, and the method comprises the following steps of: and reading the stored XY plane weight matrix, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data, determining an output image sequence after the current iteration according to the projection data, taking the output image sequence after the current iteration as a final CT reconstruction image sequence if a preset iteration stop condition is met, otherwise, taking the output image sequence after the current iteration as an input image sequence of the next iteration, and executing the next iteration, wherein the operation of reading the XY plane weight matrix is used for replacing the operation of the original generation matrix in the iteration process, so that the iteration time consumption is reduced, and the reconstruction speed is improved.
It is to be understood that both the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the disclosure.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the specification and together with the description, serve to explain the principles of the specification.
Fig. 1 is a flowchart illustrating a CT iterative reconstruction method according to an embodiment of the present invention.
Fig. 2 is a schematic diagram of memory layout rearrangement according to an embodiment of the present invention.
Fig. 3 is a functional block diagram of a CT iterative reconstruction device according to an embodiment of the present invention.
Fig. 4 is a hardware configuration diagram of an electronic device according to an embodiment of the present invention.
Detailed Description
Reference will now be made in detail to exemplary embodiments, examples of which are illustrated in the accompanying drawings. When the following description refers to the accompanying drawings, the same numbers in different drawings refer to the same or similar elements, unless otherwise indicated. The implementations described in the following exemplary examples do not represent all implementations consistent with the invention. Rather, they are merely examples of apparatus and methods consistent with aspects of embodiments of the invention as detailed in the accompanying claims.
The terminology used in the embodiments of the invention is for the purpose of describing particular embodiments of the invention only and is not intended to be limiting of embodiments of the invention. As used in this application and the appended claims, the singular forms "a," "an," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It should also be understood that the term "and/or" as used herein refers to and encompasses any or all possible combinations of one or more of the associated listed items.
It should be understood that although the terms first, second, third, etc. may be used in embodiments of the present invention to describe various information, these information should not be limited to these terms. These terms are only used to distinguish one type of information from another. For example, the first information may also be referred to as second information, and similarly, the second information may also be referred to as first information, without departing from the scope of embodiments of the present invention. The word "if" as used herein may be interpreted as "at … …" or "at … …" or "responsive to a determination", depending on the context.
Each iteration in the CT iterative reconstruction process requires recalculation of the forward projections, especially 3D (3-dimensional) forward projections, which is time consuming, resulting in a slow reconstruction speed.
The calculated 3D forward projection may be expressed mathematically as a multiplication of a three dimensional matrix, the elements of which represent the weights of the 3D forward projection, and a vector, which is the volume data required for the 3D forward projection. In the related art, in each iteration process, a three-dimensional matrix is recalculated, and then the three-dimensional matrix is multiplied by a vector of the volume data to obtain projection data.
The CT iterative reconstruction method is described in detail by way of examples.
In this embodiment, before performing iterative reconstruction, a three-dimensional weight matrix required by three-dimensional forward projection is decomposed into an XY plane weight matrix and a Z-direction weight vector, and the XY plane weight matrix is determined offline in advance and stored
Fig. 1 is a flowchart illustrating a CT iterative reconstruction method according to an embodiment of the present invention. As shown in fig. 1, in the CT iterative reconstruction method of the present embodiment, each iterative reconstruction process may include:
s101, reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data.
S102, determining an output image sequence after the iteration according to the projection data.
S103, if a preset iteration stop condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
The input image sequence in the first iteration process may be an image sequence obtained by reconstructing raw data obtained by CT scanning according to a non-iterative reconstruction algorithm.
The weights of the 3D forward projection can be decomposed into XY plane weights and Z direction weights. In the embodiment of the invention, before iterative reconstruction, an XY plane weight matrix is calculated offline in advance and stored (for example, stored in a hard disk), and in each iterative reconstruction process, the stored XY plane weight matrix is directly read, and the weight in the Z direction can be calculated in real time by utilizing CUDA (Compute Unified Device Architecture, unified computing device architecture), so that 3D forward projection calculation is completed.
Herein, use A M×N Representing an XY plane weight matrix. Offline computing A M×N The calculation time is not needed to be considered too much, so that a forward projection model calculation A with more accurate modeling can be selected M×N . The number of rows of the weight matrix is m=viewnumperrotation×channelnum, and the number of columns is n=volume_x×volume_y, because of matrix a M×N The data can be stored in a sparse matrix format to a hard disk, so that the space occupation of the hard disk can be reduced. Matrix A M×N Element a of (a) i,j Representing the weight of the jth pixel to the ith ray, an appropriate forward projection model, such as distance driving, can be selected and calculated.
In this embodiment, the distance SP (i, j) of the focal spot along the ray i to the pixel j can also be calculated off-line and stored, SP (i, j) being related to i, j only. The SP (i, j) is calculated offline and stored in the hard disk, so that the time consumption of the iterative reconstruction process can be further reduced, and the reconstruction speed is improved.
In application, an XY plane weight matrix can be calculated offline, and the calculated XY plane weight matrix is stored in a CT system. For any subject, the stored XY plane weight matrix is used as known data in the CT iterative reconstruction process, and only reading operation is needed without recalculation. In this way, the time consumption of each iteration process is reduced, so that the time consumption of the whole reconstruction process is reduced, and the reconstruction speed is improved.
It should be noted that, in a CT system, the XY plane weight matrix is the same for all subjects. Therefore, the XY plane weight matrix only needs to be calculated offline once, and the time used for the CT iterative reconstruction does not include the time for calculating the XY plane weight matrix offline.
Therefore, the embodiment avoids recalculating the 3-dimensional weight matrix in each iteration process by calculating and storing the XY plane weight matrix offline in advance, reduces the time consumption of three-dimensional forward projection in the iteration process, and improves the reconstruction speed.
In this embodiment, determining the output image sequence after the current iteration according to the projection data may include:
calculating the difference between the projection data and the corresponding CT scanning raw data to obtain error projection;
performing back projection on the error projection to obtain an error image sequence;
subtracting the corresponding image in the error image sequence from the image of the input image sequence of the current iteration to obtain an output image sequence.
In one example, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data may include:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
The formula for performing image integration on the input image sequence along the Z direction can be shown as formula (1):
for a pair of
Where intelvolome (x, y, z) represents the pixel value of the integral image, and (x, y, z) represents the z-th image, y-th row, and x-th column.
Volume (x, y, z) is the image sequence input to the 3D forward projection, i.e. the input image sequence, (x, y, z) represents the z-th image, y-th row, x-th column.
In one example, acquiring the Z-directional one-dimensional distance-driven forward projection integral corresponding to the element according to the integral image sequence may include:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
Wherein, the image integral value corresponding to the maximum value can be obtained by linear interpolation between the first image integral value and the second image integral value; the image integration value corresponding to the minimum value may be obtained by linear interpolation between the third image integration value and the fourth image integration value.
For the s-th row of detectors, the i-th ray that follows the center plane, the j-th pixel (ray i passes through pixel j), assume that the integral range of the Z-direction one-dimensional distance-driven forward projection is denoted as [ Z ] min (i,j,s),z max (i,j,s)]Wherein z is min (i, j, s) is the minimum of the integral range, z max (i, j, s) is the maximum of the integral range, then z min (i, j, s) and z max The calculation formulas of (i, j, s) are shown as formula (2) and formula (3), respectively:
z min (i,j,s)=SD/SP(i,j)*(s-SliceNum/2)*SliceThick/Interval (2)
z max (i,j,s)=SD/SP(i,j)*(s-SliceNum/2+1)*SliceThick/Interval (3)
where SD is the focal to detector distance, slicethin is the thickness of the CT detector row, sliceNum is the number of CT detector rows, and Interval is the Interval of the image sequence input to the 3D forward projection.
InteVolume (x, y, z) in z can be obtained by linear interpolation min (i, j, s) and z max Values of (i, j, s): inteVolume (x, y, z) min (i, j, s)) (i.e., the image integrated value corresponding to the minimum value) and intelvolme (x, y, z) max (i, j, s)) (i.e., the image corresponding to the maximum valueIntegral value). The Z-directed 1D distance driven forward projection integral is equal to the intelvolme (x, y, Z max (i,j,s))-InteVolume(x,y,z min (i,j,s))。
Here, it is exemplified how the image integrated value corresponding to the maximum value and the image integrated value corresponding to the minimum value are obtained by linear interpolation.
Assuming that the integral image sequence has 10 images F1-F10, and Z coordinate values corresponding to the 10 images are sequentially 1, 2, 3 and … …; z min (i,j,s)=2.5,z max (i, j, s) =7.2. For z max (i, j, s) =7.2, and Z corresponds to Z in Z coordinate values corresponding to the images F1 to F10 max (i, j, s) =7.2 adjacent values are z=7 (corresponding to image F7) and z=8 (corresponding to image F8), then element a is found i,j Image integration value V7 and element a corresponding on image F7 i,j Corresponding image integration value V8 on image F8, then z max The image integrated value corresponding to (i, j, s) can be obtained by linear interpolation between V7 and V8. Similarly, for z min (i, j, s) =2.5, and Z corresponds to Z in Z coordinate values corresponding to the images F1 to F10 min (i, j, s) =2.5 adjacent values are z=2 (corresponding to image F2) and z=3 (corresponding to image F3), then element a is found i,j Image integration value V2 and element a corresponding on image F2 i, j corresponds to the image integrated value V3 on the image F3, z min The image integrated value corresponding to (i, j, s) can be obtained by linear interpolation between V2 and V3.
Let P be v,c,s Representing the v-th projection, the c-th detector channel, the 3D forward projection values of the s-th row of detectors, P v,c,s The corresponding projection values are expressed as formula (4):
wherein ρ(s) represents the angle between the cone beam ray and the cone beam center plane, and ρ(s) is related to the detector rank s; a, a i,j For XY weight matrix A M×N Represents the weight of the jth pixel to the ith ray.
P v,c,s Corresponding projection values are accumulated to obtainTo the final projection value, the cumulative formula is shown in the following formula (5):
wherein, the relation between subscripts v and c and subscript i is: i=v. channelnum+c. Thus traversing A M×N After all elements of (3) a round of 3D forward projection calculations is completed.
Because of P v,c,s The method is positioned in the global video memory and used for storing the calculation result of the 3D forward projection, and the formula (5) needs to read and write the global video memory frequently, so that the method is time-consuming.
In one example, driving the forward projection integral according to the Z-direction one-dimensional distance corresponding to all elements in the XY plane weight matrix and the elements in the XY plane weight matrix to obtain projection data may include:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
According to the embodiment, the Z-direction 1D distance driving forward projection is subjected to block calculation, so that the cache hit rate is improved, the access to the global video memory is reduced, and the processing time is saved.
In the embodiment, a block calculation mode is adopted for the Z-direction 1D distance driving forward projection, so that the access to the global video memory is reduced. The block calculation process of this embodiment may be:
for P v,c,s And s is more than or equal to 0 and less than or equal to SliceNum-1, distributing TileSize threads for calculation, setting the x component of the Block for executing the kernel function as TileSize, and dividing the calculation into TileNum blocks.Thus, all projection components of TileNum projection values are respectively accumulated to TileNum register variables PTemp 0 ~PTemp TileNum-1 After the accumulation is finished, further adding PTemp 0 ~PTemp TileNum-1 Assignment of values to P v,c,s And a corresponding global video memory.
In one example, the sequence of integral images is stored in memory in a Z-directed preferential storage.
In this embodiment, the memory layout is rearranged, and the efficiency of CUDA combined access to the video memory can be improved through the memory layout rearrangement. Fig. 2 is a schematic diagram of memory layout rearrangement according to an embodiment of the present invention. In fig. 2, the left-hand diagram shows an integral image sequence intel Volume (X, Y, Z), which is stored in a row-first manner in the related art [ volume_x, volume_y, volume_z+1], as shown in the upper diagram on the right-hand side of fig. 2. In this embodiment, the Z-direction priority storage [ volume_z+1, volume_x, volume_y ] is changed as shown in the lower right-hand diagram of fig. 2.
According to the CT iterative reconstruction method provided by the embodiment of the invention, the three-dimensional weight matrix required by three-dimensional forward projection is regarded as the product of the XY plane weight matrix and the Z direction weight vector, the XY plane weight matrix is determined in advance in an off-line mode and stored, and each iterative reconstruction process comprises the following steps: and reading the stored XY plane weight matrix, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data, determining an output image sequence after the current iteration according to the projection data, taking the output image sequence after the current iteration as a final CT reconstruction image sequence if a preset iteration stop condition is met, otherwise, taking the output image sequence after the current iteration as an input image sequence of the next iteration, and executing the next iteration, wherein the operation of reading the XY plane weight matrix is used for replacing the operation of the original generation matrix in the iteration process, so that the iteration time consumption is reduced, and the reconstruction speed is improved.
Based on the method embodiment, the embodiment of the invention also provides a corresponding device, equipment and storage medium embodiment.
Fig. 3 is a functional block diagram of a CT iterative reconstruction device according to an embodiment of the present invention. And regarding a three-dimensional weight matrix required by three-dimensional forward projection as the product of an XY plane weight matrix and a Z-direction weight vector, and storing the XY plane weight matrix which is determined offline in advance in the device. As shown in fig. 3, in this embodiment, the CT iterative reconstruction apparatus may include:
the projection module 310 is configured to read a stored XY plane weight matrix in each iterative reconstruction process, and perform three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
a determining module 320, configured to determine an output image sequence after the current iteration according to the projection data;
the output module 330 is configured to take the output image sequence after the current iteration as a final CT reconstructed image sequence if a preset iteration stop condition is satisfied; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
In one example, the projection module 310 may be specifically configured to:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
In one example, acquiring the Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence includes:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
In one example, driving the forward projection integral according to the Z-direction one-dimensional distance corresponding to all elements in the XY plane weight matrix and the elements in the XY plane weight matrix to obtain projection data includes:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
In one example, the sequence of integral images is stored in memory in a Z-directed preferential storage.
The embodiment of the invention also provides CT equipment. Fig. 4 is a hardware configuration diagram of an electronic device according to an embodiment of the present invention. The electronic device may be, for example, a console device in a CT system. And regarding a three-dimensional weight matrix required by three-dimensional forward projection as a product of an XY plane weight matrix and a Z-direction weight vector, wherein the XY plane weight matrix which is determined in advance in an offline manner is stored in the electronic equipment. As shown in fig. 4, the electronic device includes: an internal bus 401, and a memory 402, a processor 403, and an external interface 404 connected by the internal bus, wherein:
the memory 402 is configured to store machine-readable instructions corresponding to CT iterative reconstruction logic;
the processor 403 is configured to read the machine readable instructions on the memory 402 and execute the instructions to implement the following operations:
the iterative reconstruction process comprises the following steps:
reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
determining an output image sequence after the iteration according to the projection data;
if the preset iteration stopping condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
In one example, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data, including:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
In one example, acquiring the Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence includes:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
In one example, driving the forward projection integral according to the Z-direction one-dimensional distance corresponding to all elements in the XY plane weight matrix and the elements in the XY plane weight matrix to obtain projection data includes:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
In one example, the sequence of integral images is stored in memory in a Z-directed preferential storage.
The embodiment of the invention also provides a CT system, which comprises a CT device and an electronic device, wherein:
the CT equipment is used for carrying out CT scanning on a subject to obtain raw data:
regarding a three-dimensional weight matrix required by three-dimensional forward projection as a product of an XY plane weight matrix and a Z-direction weight vector, wherein the XY plane weight matrix which is determined in advance in an offline manner is stored in electronic equipment, and the electronic equipment is used for:
the iterative reconstruction process comprises the following steps:
reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
determining an output image sequence after the iteration according to the projection data;
if the preset iteration stopping condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
In one example, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data, including:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
In one example, acquiring the Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence includes:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
In one example, driving the forward projection integral according to the Z-direction one-dimensional distance corresponding to all elements in the XY plane weight matrix and the elements in the XY plane weight matrix to obtain projection data includes:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
In one example, the sequence of integral images is stored in memory in a Z-directed preferential storage.
The embodiment of the invention also provides a computer readable storage medium, on which a computer program is stored, wherein the three-dimensional weight matrix required by three-dimensional forward projection is regarded as the product of an XY plane weight matrix and a Z direction weight vector, the computer readable storage medium also stores the XY plane weight matrix which is determined in advance in an offline manner, and the program when executed by a processor realizes the following operations:
the iterative reconstruction process comprises the following steps:
reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
determining an output image sequence after the iteration according to the projection data;
if the preset iteration stopping condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, the output image sequence after the current iteration is used as the input image sequence of the next iteration, and the next iteration is executed.
In one example, performing three-dimensional forward projection on the input image sequence by using the XY plane weight matrix to obtain projection data, including:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
In one example, acquiring the Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence includes:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
In one example, driving the forward projection integral according to the Z-direction one-dimensional distance corresponding to all elements in the XY plane weight matrix and the elements in the XY plane weight matrix to obtain projection data includes:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
In one example, the sequence of integral images is stored in memory in a Z-directed preferential storage.
For the device and apparatus embodiments, reference is made to the description of the method embodiments for the relevant points, since they essentially correspond to the method embodiments. The apparatus embodiments described above are merely illustrative, wherein the modules illustrated as separate components may or may not be physically separate, and the components shown as modules may or may not be physical, i.e., may be located in one place, or may be distributed over a plurality of network modules. Some or all of the modules may be selected according to actual needs to achieve the purposes of the present description. Those of ordinary skill in the art will understand and implement the present invention without undue burden.
The foregoing describes specific embodiments of the present disclosure. Other embodiments are within the scope of the following claims. In some cases, the actions or steps recited in the claims can be performed in a different order than in the embodiments and still achieve desirable results. In addition, the processes depicted in the accompanying figures do not necessarily require the particular order shown, or sequential order, to achieve desirable results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.
Other embodiments of the present description will be apparent to those skilled in the art from consideration of the specification and practice of the invention disclosed herein. This specification is intended to cover any variations, uses, or adaptations of the specification following, in general, the principles of the specification and including such departures from the present disclosure as come within known or customary practice within the art to which the specification pertains. It is intended that the specification and examples be considered as exemplary only, with a true scope and spirit of the specification being indicated by the following claims.
It is to be understood that the present description is not limited to the precise arrangements and instrumentalities shown in the drawings, which have been described above, and that various modifications and changes may be made without departing from the scope thereof. The scope of the present description is limited only by the appended claims.
The foregoing description of the preferred embodiments is provided for the purpose of illustration only, and is not intended to limit the scope of the disclosure, since any modifications, equivalents, improvements, etc. that fall within the spirit and principles of the disclosure are intended to be included within the scope of the disclosure.

Claims (8)

1. A CT iterative reconstruction method is characterized in that a three-dimensional weight matrix required by three-dimensional forward projection is regarded as the product of an XY plane weight matrix and a Z direction weight vector, and the XY plane weight matrix is determined in advance in an off-line mode and stored;
the iterative reconstruction process comprises the following steps:
reading a stored XY plane weight matrix, and performing three-dimensional forward projection on an input image sequence by using the XY plane weight matrix to obtain projection data;
determining an output image sequence after the iteration according to the projection data;
if the preset iteration stopping condition is met, taking the output image sequence after the iteration as a final CT reconstructed image sequence; otherwise, taking the output image sequence after the current iteration as the input image sequence of the next iteration, and executing the next iteration;
the three-dimensional forward projection is performed on the input image sequence by using the XY plane weight matrix to obtain projection data, and the method comprises the following steps:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
2. The method of claim 1, wherein obtaining a Z-directed one-dimensional distance-driven forward projection integral corresponding to the element from the sequence of integral images comprises:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
3. The method of claim 1, wherein driving forward projection integrals according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and elements in the XY plane weight matrix to obtain projection data comprises:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
4. The method of claim 1, wherein the sequence of integrated images is stored in memory in a Z-directed preferential storage.
5. A CT iterative reconstruction apparatus, wherein a three-dimensional weight matrix required for three-dimensional forward projection is regarded as a product of an XY plane weight matrix and a Z-direction weight vector, the apparatus storing therein the XY plane weight matrix determined offline in advance, the apparatus comprising:
the projection module is used for reading the stored XY plane weight matrix in each iterative reconstruction process, and carrying out three-dimensional forward projection on the input image sequence by utilizing the XY plane weight matrix to obtain projection data;
the determining module is used for determining an output image sequence after the iteration according to the projection data;
the output module is used for taking the output image sequence after the iteration as a final CT reconstructed image sequence if a preset iteration stop condition is met; otherwise, taking the output image sequence after the current iteration as the input image sequence of the next iteration, and executing the next iteration;
the projection module is specifically used for:
performing image integration on the input image sequence along the Z direction to obtain an integrated image sequence;
traversing each detector row for each element of the XY plane weight matrix, and acquiring a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the integral image sequence;
and driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and obtaining projection data according to the elements in the XY plane weight matrix.
6. The apparatus of claim 5, wherein obtaining a Z-directed one-dimensional distance-driven forward projection integral corresponding to the element from the sequence of integral images comprises:
for each detector row s, reading the distance SP (i, j) from ray i passing through pixel j to pixel j; wherein the distance SP (i, j) is obtained and stored by offline calculation in advance;
determining the maximum value and the minimum value of the integral range of the Z-direction one-dimensional distance driving forward projection according to preset detector parameters and the distance SP (i, j);
determining a first value and a second value adjacent to the maximum value from the Z coordinate values corresponding to the integrated image sequence, and determining a third value and a fourth value adjacent to the minimum value;
obtaining an image integral value corresponding to the maximum value according to a first image integral value corresponding to the element on the integral image corresponding to the first value and a second image integral value corresponding to the element on the integral image corresponding to the second value; obtaining an image integral value corresponding to the minimum value according to a third image integral value corresponding to the element on the integral image corresponding to the third value and a fourth image integral value corresponding to the element on the integral image corresponding to the fourth value;
and determining a Z-direction one-dimensional distance driving forward projection integral corresponding to the element according to the image integral corresponding to the maximum value and the image integral corresponding to the minimum value.
7. The apparatus of claim 5, wherein driving forward projection integration according to the Z-direction one-dimensional distances corresponding to all elements in the XY plane weight matrix and elements in the XY plane weight matrix to obtain projection data comprises:
distributing the calculation of the projection data to TileSize threads, and processing projection data calculation corresponding to TileNum detector rows by each thread; the TileNum detectors corresponding to each thread belong to the TileNum groups;
respectively accumulating the calculation results of the TileSize threads into TileNum register variables according to the corresponding groups;
after the accumulation is finished, the accumulation results of the TileNum register variables are written into the global video memory corresponding to the projection data p (v, c, s).
8. The apparatus of claim 5, wherein the sequence of integrated images is stored in memory in a Z-directed preferential storage.
CN202110443876.3A 2021-04-23 2021-04-23 CT iterative reconstruction method and device Active CN113284205B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110443876.3A CN113284205B (en) 2021-04-23 2021-04-23 CT iterative reconstruction method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110443876.3A CN113284205B (en) 2021-04-23 2021-04-23 CT iterative reconstruction method and device

Publications (2)

Publication Number Publication Date
CN113284205A CN113284205A (en) 2021-08-20
CN113284205B true CN113284205B (en) 2024-01-02

Family

ID=77277262

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110443876.3A Active CN113284205B (en) 2021-04-23 2021-04-23 CT iterative reconstruction method and device

Country Status (1)

Country Link
CN (1) CN113284205B (en)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105931280A (en) * 2016-03-29 2016-09-07 中北大学 GPU-based fast three-dimensional CT (computed tomography) iterative reconstruction method
CN111110260A (en) * 2019-12-24 2020-05-08 沈阳先进医疗设备技术孵化中心有限公司 Image reconstruction method and device and terminal equipment
CN111192228A (en) * 2020-01-02 2020-05-22 沈阳先进医疗设备技术孵化中心有限公司 Image processing method and device, CT (computed tomography) equipment and CT system

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7254209B2 (en) * 2003-11-17 2007-08-07 General Electric Company Iterative CT reconstruction method using multi-modal edge information
US9508163B2 (en) * 2013-06-14 2016-11-29 General Electric Company Accelerated iterative reconstruction
EP3014577A1 (en) * 2013-06-28 2016-05-04 Koninklijke Philips N.V. Methods for generation of edge-preserving synthetic mammograms from tomosynthesis data
JP6215449B2 (en) * 2014-03-14 2017-10-18 株式会社日立製作所 X-ray CT apparatus and processing apparatus
JP6766045B2 (en) * 2014-11-20 2020-10-07 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. How to generate a synthetic mammogram from tomosynthesis data
US9911208B2 (en) * 2016-04-11 2018-03-06 Toshiba Medical Systems Corporation Apparatus and method of iterative image reconstruction using regularization-parameter control

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105931280A (en) * 2016-03-29 2016-09-07 中北大学 GPU-based fast three-dimensional CT (computed tomography) iterative reconstruction method
CN111110260A (en) * 2019-12-24 2020-05-08 沈阳先进医疗设备技术孵化中心有限公司 Image reconstruction method and device and terminal equipment
CN111192228A (en) * 2020-01-02 2020-05-22 沈阳先进医疗设备技术孵化中心有限公司 Image processing method and device, CT (computed tomography) equipment and CT system

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
J. J. Scheins et al.Fully-3D PET Image Reconstruction Using Scanner-Independent, Adaptive Projection Data and Highly Rotation-Symmetric Voxel Assemblies.《IEEE Transactions on Medical Imaging 》.2011,第30卷(第3期),第879 - 892页. *
赵星 ; 胡晶晶 ; 张慧滔 ; 张朋 ; .基于图形处理器的锥束CT快速迭代重建算法.北京理工大学学报.2010,第30卷(第11期),第1310-1314页. *

Also Published As

Publication number Publication date
CN113284205A (en) 2021-08-20

Similar Documents

Publication Publication Date Title
JP3527796B2 (en) High-speed three-dimensional image generating apparatus and method
CN108896943A (en) A kind of magnetic resonance quantitative imaging method and device
US8983162B2 (en) Method and apparatus for estimating monte-carlo simulation gamma-ray scattering in positron emission tomography using graphics processing unit
US8620054B2 (en) Image reconstruction based on accelerated method using polar symmetries
CN109377500B (en) Image segmentation method based on neural network and terminal equipment
US9189870B2 (en) Method, computer readable medium and system for tomographic reconstruction
Liu et al. GPU-based branchless distance-driven projection and backprojection
CN110811667B (en) High-precision PET reconstruction method and device based on GPU acceleration
Flores et al. CT image reconstruction based on GPUs
Cuomo et al. A GPU parallel implementation of the local principal component analysis overcomplete method for DW image denoising
KR101283266B1 (en) Method and apparatus for monte-carlo simulation gamma-ray scattering estimation in positron emission tomography using graphics processing unit
Corda-D’Incan et al. Memory-efficient training for fully unrolled deep learned PET image reconstruction with iteration-dependent targets
CN113284205B (en) CT iterative reconstruction method and device
US8416239B2 (en) Intermediate image generation method, apparatus, and program
CN106910157A (en) The image rebuilding method and device of a kind of multistage parallel
DE102022120819A1 (en) QUANTIZED NEURAL NETWORK TRAINING AND INFERENCE
Zhao et al. Fast projection algorithm for voxel arrays with object dependent boundaries
CN115375583A (en) PET parameter image enhancement method, device, equipment and storage medium
CN107730464A (en) Image noise reduction parallel algorithm based on Block- matching
Sampson et al. Investigating multi-threaded SIMD for helical CT reconstruction on a CPU
US6992704B2 (en) Image processing apparatus and method, recording medium, and program
CN116739933A (en) PET image generation method integrating system response matrix and neural network
CN111161366A (en) Image reconstruction method and device, terminal equipment and storage medium
CN112991482B (en) GPU-based rapid reconstruction imaging method and device and readable storage medium
CN116091219A (en) Single-obstacle European option Vega numerical calculation method and device

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20240204

Address after: 110167 No. 177-1 Innovation Road, Hunnan District, Shenyang City, Liaoning Province

Patentee after: Shenyang Neusoft Medical Systems Co.,Ltd.

Country or region after: China

Address before: Room 336, 177-1, Chuangxin Road, Hunnan New District, Shenyang City, Liaoning Province

Patentee before: Shenyang advanced medical equipment Technology Incubation Center Co.,Ltd.

Country or region before: China