BIOMEDICAL DATA ANALYSIS ON HETEROGENEOUS PLATFORM
Methods and apparatus for biomedical data analysis to produce enhanced images of tubular structures are disclosed. A Gaussian convolution of an input image is used to calculate a Hessian matrix. An Eigen decomposition of the Hessian matrix produces eigenvectors and eigenvalues, which are sorted to determine bright tubular structure detection according to high and low values that represent brightness, and structure shape. A tubularity computation calculates the probability of a voxel of interest being part of a tubular network. Embodiments may be implemented to share computer resources such as between a computer processing unit (CPU) and a graphic processing unit (GPU).
This application claims the benefit of U.S. Provisional Patent Application No. 61/657,451 filed Jun. 8, 2012, which is incorporated by reference as if fully set forth.
FIELD OF THE INVENTIONThe present invention is generally directed to graphical enhancement of images including, but not limited to, biomedical images such as those related to vascular bodies and vasculature extraction.
BACKGROUNDBiomedical image generation and analysis involves detection and extraction of vascular structures and tissues from a two dimensional images or three dimensional renderings. Computations required to detect and extract such structures are intensive and can have undesirable latencies. Known techniques that can produce images with adequate fine detail to properly extract fine structures such as vascular structures are not capable of real time computations, or suffer significant latency.
SUMMARY OF EMBODIMENTSMethods and apparatus for biomedical data analysis to produce enhanced images of tubular or tubularlike structures, such as vascular bodies, are disclosed. A Gaussian convolution of an input image is used to calculate a Hessian matrix. An Eigen decomposition of the Hessian matrix produces eigenvectors and eigenvalues, which are sorted to determine bright tubular structure detection according to high and low values that represent brightness, and structure shape. A tubularity computation, also referred to as a vesselness computation, calculates the probability of a voxel of interest being part of a tubular network, such as a vascular network. Embodiments may be implemented to share computer resources such as between a computer processing unit (CPU) and a graphic processing unit (GPU).
In one embodiment, a method for image enhancement is provided. A Gaussian convolution of voxels of an input image is performed. A Hessian matrix of the convolution is calculated. An Eigen decomposition of the Hessian matrix is performed to derive eigenvectors and eigenvalues of the voxels. The eigenvectors are sorted according to the magnitude of the corresponding eigenvalues. A tubularity computation, also referred to as a vesselness computation, is performed to calculate the probability of a voxel of the image being part of a tubular structure representation in the image whereby voxels determined to be part of tubular structures, such as vascular structures, are used to provide an enhanced depiction of the structures.
Such an embodiment may be implemented in stages where the Gaussian convolution is performed a first processing stage, the a Hessian matrix is calculated in a second processing stage, the Eigen decomposition and sorting the eigenvectors is performed in a third processing stage, and the a vesselness computation is performed in a fourth processing stage. A plurality of computer resources of a device may share in the processing of the voxel data. For example, one computer resource, such as a central processing unit (CPU), may perform overall tracking and coordination of the results between stages for voxels of the image and another computer resource, such as a graphics processing unit (GPU), may perform one or more of the processing stages.
In another embodiment an apparatus for biomedical image enhancement is provided. A first computational component is configured to perform a Gaussian convolution of voxels of an input image. A second computational component is configured to calculate a Hessian matrix of the convolution. A third computational component is configured to perform an Eigen decomposition of the Hessian matrix to derive eigenvectors and eigenvalues of the voxels and to sort the eigenvectors according to the magnitude of the corresponding eigenvalues. A fourth computational is component configured to perform a tubularity computation, also referred to as a vesselness computation, to calculate the probability of a voxel of the image being part of a tubular structure representation in the image whereby voxels determined to be part of tubular structures, such as vascular structures, are used to provide an enhanced depiction of the structures.
In one example, the four computational components are implemented in a central processing unit (CPU) of a hand held battery powered device.
Embodiments include devices having a plurality of computer resources wherein one computer resource, such as a central processing unit (CPU), is configured to perform overall tracking and coordination of the results between the four computational components and one or more of the four computational components are implemented in a different computer resource, such as a graphics processing unit (GPU). Embodiments also include devices having a central processing unit (CPU) configured to perform overall tracking and coordination of the results between the four computational components and further configured to selectively implement the four computational components.
Embodiments may also include an input configured to receive the voxels of the input image and an output configured to output voxels of the processed image.
In further embodiments, computer readable medium having stored thereon a program that when executed perform steps of the method embodiments are provided. The computer readable medium may implement a program where the Gaussian convolution is performed in a first processing stage, the Hessian matrix is calculated in a second processing stage, the Eigen decomposition and sorting the eigenvectors is performed in a third processing stage, and the vesselness computation is performed in a fourth processing stage wherein the program controls a central processing unit (CPU) to perform overall tracking and coordination of the results between stages for voxels of the image and a graphics processing unit (GPU) to perform at least one of the processing stages.
A more detailed understanding may be had from the following description, given by way of example in conjunction with the accompanying drawings wherein:
Methods and apparatus are used to perform image data analysis according to the embodiments described herein. From the data analysis, images produced by, for example, scanned organs and tissues are analyzed for detection and extraction of tubular structures, such as blood vessels.
The processor 102 may include a central processing unit (CPU), a graphics processing unit (GPU), a CPU and GPU located on the same die, or one or more processor cores, wherein each processor core may be a CPU or a GPU. The memory 104 may be located on the same die as the processor 102, or may be located separately from the processor 102. The memory 104 may include a volatile or nonvolatile memory, for example, random access memory (RAM), dynamic RAM, or a cache. The processor 102 may execute the biomedical data analysis method which may be stored as executable code on the memory 104.
The storage 106 may include a fixed or removable storage, for example, a hard disk drive, a solid state drive, an optical disk, or a flash drive. The input devices 108 may include a keyboard, a keypad, a touch screen, a touch pad, a detector, a microphone, an accelerometer, a gyroscope, a biometric scanner, or a network connection (e.g., a wireless local area network card for transmission and/or reception of wireless IEEE 802 signals). The output devices 110 may include a display, a speaker, a printer, a haptic feedback device, one or more lights, an antenna, or a network connection (e.g., a wireless local area network card for transmission and/or reception of wireless IEEE 802 signals).
The input driver 112 communicates with the processor 102 and the input devices 108, and permits the processor 102 to receive input from the input devices 108. The output driver 114 communicates with the processor 102 and the output devices 110, and permits the processor 102 to send output to the output devices 110. It is noted that the input driver 112 and the output driver 114 are optional components, and that the device 100 will operate in the same manner if the input driver 112 and the output driver 114 are not present.
In a first stage 210, a Gaussian convolution is performed on an input image I, using three onedimensional (1D) kernels instead of one threedimensional (3D) kernel. This greatly reduces the overall complexity of the image processing. From the convolution, a Hessian matrix H is computed in a second stage 220. For the third stage 230, an Eigen decomposition of the Hessian matrix computed in the second stage 220 is performed to derive its eigenvectors and eigenvalues, followed by sorting the eigenvectors according to the magnitude of the corresponding eigenvalues. Finally, in a fourth stage 240, a tubularity computation, also referred to as a vesselness computation, calculates the probability of a voxel being part of a tubular network in order to provide the enhanced detail of a tubular structure, such as a vessel structure, in the processed image. The tubularity value, i.e. vesselness value, of each voxel is then used as intensity value of each corresponding voxel in the enhanced image, only with linear scaling if necessary, to produce the postprocessed image, for example,
The Gaussian convolution performed in the first stage 210 may be represented by the equation:
L=((I*G_{x}(σ))*G_{y}(σ))*G_{z}(σ) Equation (1)
where I is the input image, G_{x}, G_{y}, G_{z }are 1D Gaussian kernels in 3D space in dimensions x, y and z, and σ is the size of the Gaussian kernel. The Hessian matrix computation of the convolution L performed in the second stage 220 may be represented as:
For the third stage 230, an Eigen decomposition of such computed Hessian matrix is performed to derive its eigenvectors and eigenvalues. The subsequent sorting of the eigenvectors according to the magnitude of the corresponding eigenvalues permits analysis of the main modes of secondorder variation in image intensity to determine the type of local structure in different portions of the image.
The diagram shown in
Where gaps in the estimation of the tubular structure segment (e.g., a vascular segment) occur due to errors in the calculation, tensor voting based on mathematical diffusion may be applied. The longest connected components, such as the rendered vessels illustrated in
Each stage of the method may utilize multithreaded or single threaded execution on a CPU, a GPU or a combination of both CPU and GPU. For example,
With reference to the steps shown in
For this example, the shortest observed execution time was associated with GPU execution of all four computational stages in connection with overall CPU control of the processing. However, computational times will vary with the extent of the tubular structures within the subject image. The four computational stages as shown in
The processing of the 3D Gaussian convolution of 3D images in 1D kernels consecutively for each dimension, instead of using one 3D kernel, greatly reduces the computational complexity of the overall processing. This is applicable to other numbers of dimensionalities. Coupling the reduced processing complexity with selective work sharing between computing resources such as CPUs and GPUs makes implementation of the method practical for smaller and handheld devices since the image enhancement is both fast and requires less power than attempting to provide such enhanced images through a more brute force computational process.
With reference to
ℑ(x)=({right arrow over (υ)}_{2}·∇I)^{2}+({right arrow over (υ)}_{3}·∇I)^{2}≈0
with ridge points of a vessel satisfying the following:
{right arrow over (υ)}_{2}·∇I=0 {right arrow over (υ)}_{3}·∇I=0 λ_{2}<0 λ_{3}<0
This includes a ridge detection and traversal calculation illustrated in
A tubular extraction method, such as the example vasculature extraction algorithm, may apply a templatebased approach such as represented by:
An artificial tubular template is represented according to an image contrast and intensity scale, such as is exemplified in
Template fitting may be based on voxel x, center of template x_{0}, vector v for the tubular structure direction, and background intensity m, where a weighted W difference between template intensity and image intensity is determined. The value ε(x) represents an error value. This may be represented as:
Graphic examples are provided as
An example pseudo code methodology for determining vascular trees of an image S is as follows given Input as a vasculature image S and output are the vascular trees in that image.

 Step 1: Parallel for each voxel v in the image S, compute the probability p(v) representing how likely a voxel v is part of the tubular structure.
 Step 2: Sort all voxels based on p(v) value.
 Step 3: Parallel for each voxel with a high probability (above a predefined threshold) of being in the tubular tree:
 A: Extract the vascular branch from that voxel, noted as parent voxel, by processing a region around that voxel; if branching/divergence happens, parallel for each such children voxel.
 B: Extract the tubular branch (e.g. vascular branch) from each child voxel also by processing a region around the voxel.
It should be understood that many variations are possible based on the disclosure herein. Although features and elements are described above in particular combinations, each feature or element may be used alone without the other features and elements or in various combinations with or without other features and elements.
The methods provided may be implemented in a general purpose computer, a processor, or a processor core. Suitable processors include, by way of example, a general purpose processor, a special purpose processor, a conventional processor, a digital signal processor (DSP), a plurality of microprocessors, one or more microprocessors in association with a DSP core, a controller, a microcontroller, Application Specific Integrated Circuits (ASICs), Field Programmable Gate Arrays (FPGAs) circuits, any other type of integrated circuit (IC), and/or a state machine. Such processors may be manufactured by configuring a manufacturing process using the results of processed hardware description language (HDL) instructions and other intermediary data including netlists (such instructions capable of being stored on a computer readable media). The results of such processing may be maskworks that are then used in a semiconductor manufacturing process to manufacture a processor which implements aspects of the present invention.
The methods or flow charts provided herein may be implemented in a computer program, software, or firmware incorporated in a computerreadable storage medium for execution by a general purpose computer or a processor. Examples of computerreadable storage mediums include a read only memory (ROM), a random access memory (RAM), a register, cache memory, semiconductor memory devices, magnetic media such as internal hard disks and removable disks, magnetooptical media, and optical media such as CDROM disks, and digital versatile disks (DVDs).
Claims
1. A method for image enhancement comprising:
 performing a Gaussian convolution of voxels of an input image;
 calculating a Hessian matrix of the convolution;
 performing an Eigen decomposition of the Hessian matrix to derive eigenvectors and eigenvalues of the voxels;
 sorting the eigenvectors according to the magnitude of the corresponding eigenvalues; and
 performing a tubularity computation to calculate the probability of a voxel of the image being part of a tubular structure representation in the image whereby voxels determined to be part of tubular structures are used to provide an enhanced depiction of the image.
2. The method of claim 1 wherein the Gaussian convolution is performed with respect to three onedimensional kernels to produce a convolution L based on where I is the input image, Gx, Gy, Gz are onedimensional Gaussian kernels with respect to three dimensional space having dimensions x, y and z, and σ is the size of the Gaussian kernel.
 L=((I*Gx(σ))*Gy(σ))*Gz(σ)
3. The method of claim 2 wherein the Hessian matrix is calculated as a matrix H in accordance with: H = ( ∂ 2 L ∂ x ∂ x ∂ 2 L ∂ x ∂ y ∂ 2 L ∂ x ∂ z ∂ 2 L ∂ y ∂ x ∂ 2 L ∂ y ∂ y ∂ 2 L ∂ y ∂ z ∂ 2 L ∂ z ∂ x ∂ 2 L ∂ z ∂ y ∂ 2 L ∂ z ∂ z ).
4. The method of claim 2 where the Gaussian convolution is performed a first processing stage, the a Hessian matrix is calculated in a second processing stage, the Eigen decomposition and sorting the eigenvectors is performed in a third processing stage, and the tubularity computation is performed in a fourth processing stage.
5. The method of claim 4 performed by a plurality of computer resources of a device wherein one computer resource performs overall tracking and coordination of the results between stages for voxels of the image and another computer resource performs at least one of the processing stages.
6. The method of claim 4 wherein a central processing unit (CPU) performs overall tracking and coordination of the results between stages for voxels of the image and a graphics processing unit (GPU) performs at least one of the processing stages.
7. The method of claim 4 wherein a central processing unit (CPU) performs overall tracking and coordination of the results between stages for voxels of the image and a graphics processing unit (GPU) performs the four processing stages.
8. The method of claim 1 where the image is a biomedical image wherein performing a tubularity computation is performed by performing a vesselness computation to calculate the probability of a voxel of the image being part of a vascular structure representation in the image whereby voxels determined to be part of vascular structures are used to provide an enhanced depiction of the biomedical image.
9. The method of claim 8 wherein:
 the Gaussian convolution is performed with respect to three onedimensional kernels to produce a convolution L based on L=((I*Gx(σ))*Gy(σ))*Gz(σ)
 where I is the input image, Gx, Gy, Gz are onedimensional Gaussian kernels with respect to three dimensional space having dimensions x, y and z, and σ is the size of the Gaussian kernel.
 the Gaussian convolution is performed a first processing stage, the a Hessian matrix is calculated in a second processing stage, the Eigen decomposition and sorting the eigenvectors is performed in a third processing stage, and the vesselness computation is performed in a fourth processing stage.
10. The method of claim 9 performed by a plurality of computer resources of a device wherein a central processing unit (CPU) performs overall tracking and coordination of the results between stages for voxels of the image and a graphics processing unit (GPU) performs at least one of the processing stages.
11. An apparatus for image enhancement comprising:
 a first computational component configured to perform a Gaussian convolution of voxels of an input image;
 a second computational component configured to calculate a Hessian matrix of the convolution;
 a third computational component configured to perform an Eigen decomposition of the Hessian matrix to derive eigenvectors and eigenvalues of the voxels and to sort the eigenvectors according to the magnitude of the corresponding eigenvalues; and
 a fourth computational component configured to perform a tubularity computation to calculate the probability of a voxel of the image being part of a tubular structure representation in the image whereby voxels determined to be part of tubular structures are used to provide an enhanced depiction of the image.
12. The apparatus of claim 11 wherein the first computational component is configured to perform the Gaussian convolution with respect to three onedimensional kernels to produce a convolution L based on where I is the input image, Gx, Gy, Gz are onedimensional Gaussian kernels with respect to three dimensional space having dimensions x, y and z, and σ is the size of the Gaussian kernel.
 L=((I*Gx(σ))*Gy(σ))*Gz(σ)
13. The apparatus of claim 12 wherein the second computational component is configured to calculate the Hessian matrix as a matrix H in accordance with: H = ( ∂ 2 L ∂ x ∂ x ∂ 2 L ∂ x ∂ y ∂ 2 L ∂ x ∂ z ∂ 2 L ∂ y ∂ x ∂ 2 L ∂ y ∂ y ∂ 2 L ∂ y ∂ z ∂ 2 L ∂ z ∂ x ∂ 2 L ∂ z ∂ y ∂ 2 L ∂ z ∂ z ).
14. The apparatus of claim 12 wherein the four computational components are implemented in a central processing unit (CPU) of the apparatus.
15. The apparatus of claim 12 having a plurality of computer resources wherein one computer resource is configured to perform overall tracking and coordination of the results between the four computational components and at least one of the four computational components is implemented in a different computer resource.
16. The apparatus of claim 12 having a central processing unit (CPU) configured to perform overall tracking and coordination of the results between the four computational components and at least one of the four computational components in implemented in a graphics processing unit (GPU).
17. The apparatus of claim 12 having a central processing unit (CPU) configured to perform overall tracking and coordination of the results between the four computational components and the four computational components are implemented in a graphics processing unit (GPU).
18. The apparatus of claim 12 having a central processing unit (CPU) configured to perform overall tracking and coordination of the results between the four computational components and to selectively implement the four computational components.
19. The apparatus of claim 18 further comprising an input configured to receive the voxels of the input image and an output configured to output voxels of the processed image.
20. A computer readable medium having stored thereon a program that when executed performs the following method steps:
 performing a Gaussian convolution of voxels of an input image;
 calculating a Hessian matrix of the convolution;
 performing an Eigen decomposition of the Hessian matrix to derive eigenvectors and eigenvalues of the voxels;
 sorting the eigenvectors according to the magnitude of the corresponding eigenvalues; and
 performing a tubularity computation to calculate the probability of a voxel of the image being part of a tubular structure representation in the image whereby voxels determined to be part of tubular structures are used to provide an enhanced depiction of the image.
21. The computer readable medium of claim 20 wherein the Gaussian convolution is performed with respect to three onedimensional kernels to produce a convolution L based on where I is the input image, Gx, Gy, Gz are onedimensional Gaussian kernels with respect to three dimensional space having dimensions x, y and z, and σ is the size of the Gaussian kernel.
 L=((I*Gx(σ))*Gy(σ))*Gz(σ)
22. The computer readable medium of claim 21 wherein the Hessian matrix is calculated as a matrix H in accordance with: H = ( ∂ 2 L ∂ x ∂ x ∂ 2 L ∂ x ∂ y ∂ 2 L ∂ x ∂ z ∂ 2 L ∂ y ∂ x ∂ 2 L ∂ y ∂ y ∂ 2 L ∂ y ∂ z ∂ 2 L ∂ z ∂ x ∂ 2 L ∂ z ∂ y ∂ 2 L ∂ z ∂ z ).
23. The computer readable medium of claim 20 where the Gaussian convolution is performed in a first processing stage, the Hessian matrix is calculated in a second processing stage, the Eigen decomposition and sorting the eigenvectors is performed in a third processing stage, and the tubularity computation is performed in a fourth processing stage wherein the program controls a central processing unit (CPU) to perform overall tracking and coordination of the results between stages for voxels of the image and a graphics processing unit (GPU) to perform at least one of the processing stages.
24. An enhanced image produced in accordance with the method of claim 1.
25. An enhanced biomedical image produced in accordance with the method of claim 8.
26. An enhanced biomedical image produced in accordance with the method of claim 9.
Type: Application
Filed: Jun 4, 2013
Publication Date: Dec 12, 2013
Patent Grant number: 9240034
Inventor: Dongping Zhang (Austin, TX)
Application Number: 13/909,541
International Classification: G06T 5/00 (20060101);