530 likes | 708 Views
GPGPU in Medical Imaging Applications. Kevin Gorczowski. Overview. Random Walks for Interactive Organ Segmentation in Two and Three Dimensions Interactive, GPU-Based Level Sets for 3D Brain Tumor Segmentation Image Registration by a Regularized Gradient Flow
E N D
GPGPU in Medical Imaging Applications Kevin Gorczowski
Overview • Random Walks for Interactive Organ Segmentation in Two and Three Dimensions • Interactive, GPU-Based Level Sets for 3D Brain Tumor Segmentation • Image Registration by a Regularized Gradient Flow • Accelerating Popular Tomographic Reconstruction Algorithms on Commodity PC Graphics Hardware
GPUs and Medical Imaging • Computations involving 2D and 3D images • Already a natural fit for use of GPUs • Often involve numerical computations • ODEs and PDEs • Numerical accuracy not always top priority • Visually correct and speed most important
Image Segmentation • Find the boundary or interior of a specific organ or structure • Hand-tracing still a standard clinical practice • Example use: radiation cancer treatment • Take CT scan on initial day • Segment organ • Use segmentation to aim radiation • Problem: organs move between days! • Hand-tracing much too time consuming to do each day
Random Walks for Interactive Organ Segmentation in Two and Three Dimensions (2005) • Goal: Segmentation using only small amount of user interaction • User specifies seed points • Seeds labeled according to organ/structure • Segmentation propagates toward seed points
Seed and Segmentation Examples Seeds: gray Segmentations: black
Methodology • Unlabeled pixels (non-seeds) • Send out random walker • Determine probability that random walker will reach a labeled pixel first, ahead of all other labels • Assign pixel maximum probability label
Formulating Images as Weighted Graphs • Nodes: pixels • Edges: between neighboring pixels • Edge weight: function of intensity difference between pixel and its neighbor • Weight function • g: pixel intensity, β: free parameter • Able to use same β throughout
Solution • Probability of random walker reaching seed point same as solution to Dirichlet problem • Partial differential equation on interior of region with conditions at boundary of region • Boundary conditions at seed pixels • 1 if seed pixel belongs to label being searched for • 0 otherwise • Computation of solution involves solving linear system with graph Laplacian matrix
GPU Implementation • Linear system LX = -BM must be solved for each seed label • Only M, seed boundary conditions, changes • Can use RGBA channels to solve for 4 labels simultaneously • Progress of segmentation can be updated on screen
GPU Implementation • L is symmetric and diagonal is determined by off-diagonal elements • Can store L as size n = # of pixels, instead of n x n • Textures are same size as original images • Linear system solved using conjugate gradient • Only matrix-vector product and vector inner-product required
Interactive, GPU-Based Level Sets for 3D Brain Tumor Segmentation (2003) • Uses deformable model approach • Before, started with pixel seeds and labeled pixels individually • Here, start with template model of organ/structure • Deform model to fit target image • Resulting model represents boundary
Deformable Models • Two classes of deformable shape models • Parametric • Parameterized curves or surfaces • Spherical harmonics, wavelets • Geometric • Curves as embedded, implicit level sets of higher-dimensional functions • Can change topology (may or may not be advantage depending on application)
Level Sets • Curve or surface: all points such that some function φ(x) = 0 • φ: R3 -> R, x: position • Can describe motion (deformation) as PDE • v(t): pixel-wise velocities • Implement deformations by choosing v’s
Velocities for Segmentation • Velocities usually have two components when used for segmentation • Data term • Smoothness term • Data term drives curve toward boundaries • Typically areas of high contrast in image • Smoothness term constrains curve from becoming too irregular • Prevents “leaking” out of small discontinuities in image edges
Problems • Data term usually introduces free parameters • Their data term: • I: image intensity • ε: determines range of values considered inside object • T: determines how bright object is • Basically, T is mean intensity and ε is variance • User has to choose these free parameters • If segmentation runs at interactive speed (GPU), user can be much more productive
GPU Implementation Issues • For stability, curve must move at most one pixel at each iteration • Results in many iterations before solution • Need to keep as much on GPU as possible • Major speedup if only places where φ is close to zero are considered • How to pack for GPU • Changes after each iteration
Data Packing • Subdivide data into 16x16 tiles • Only tiles with non-zero derivatives sent to GPU
Data Packing • Evaluation of PDE requires discrete derivatives • Pixels need to access neighbors • CPU calculates and sends texture coordinates of neighbor pixels in packed format • GPU does neighbor lookups and finite differences used for gradient and curvature • CPU also sends vertices of active tiles for quad rendering
Data Packing • CPU needs to active tiles for each iteration so it can calculate neighbor pixels • GPU writes out a small (<64KB), encoded texture telling which tiles are active • Checks if a tile boundary has been crossed • Limits CPU <-> GPU bus use
Performance • CPU: 1.7GHz Xeon • 7-8 iterations/sec • “Highly-optimized, sparse-field, CPU-based solver" • Radeon 9700 Pro • 60-70 iteration/sec • Running time of GPU implementation scales linearly with number of active voxels • Overhead for feedback image calculation about 15% of total GPU time
Results Semi-automatic result has less discontinuities across slices than hand-traced segmentation
Image Registration by a Regularized Gradient Flow (2004) • Image registration • Try to get intensities in multiple images to match spatially • Simple case: use similarity transforms to align images as best as possible • Complex case: non-rigid registration • Each pixel has its own displacement
Application: Atlas Formation • Lorenzen (UNC) created sharp “mean” images by finding mean deformation
Registration Formulation • Must define a measure of how closely two images match • Can formulate as an energy function
Registration as Optimization • Given formulation as energy function • Problem becomes a global optimization of an objective function • Find minimum of energy function
Accelerating Popular Tomographic Reconstruction Algorithms on Commodity Graphics Hardware • Computed Tomography • CT scan, used to be CAT scan • X-ray source • Object attenuates x-rays • Collector (2D sheet) measures left-over radiation • Source and collector rotate around object
Reconstruction • Only information acquired through scan is 2D “image” at collector • Have collector image for each angle φ(position of source/collector on circle) • Must solve for attenuation of scanned object at points on 3D grid
Formulation • Amount of radiation collected at pixel (u,v) of collector for angle φ • μ: attenuation (unknown) • Q0: original x-ray energy at source • L: distance between source and collector • Integrate attenuation along x-ray
Formulation • Rewrite • i: pixel index of collector • Voxel form (loop through object voxel Grid)
Formulation • wij : weight of voxel j’s contribution to detector pixel i • Determined ahead of time by interpolation and integration rules
Solving for Attenuation • Mφ : number of pixels in collector for angle φ • Know qi and wij, solve for attenuation (backprojection)
Final Image Reconstruction • Feldkamp algorithm • SART iterative algorithm • OS-EM algorithm
GPU Implementation • Each iteration requires at least one backprojection and projection • Each backprojection is O(n3) • No real way around complexity • Must use brute force speed • Many CT machines use FPGAs or ASICs • Expensive • Inflexible
GPU Implementation • Volume data stored at stacks of textures • Discrete form of projection and backprojection operations • Updates of current volume attenuation performed through texture blends
GPU Implementation • Can use RGBA channels to compute orthogonal projections • Projection matrices are equal • Four 90° increments of φ • Can’t do this in SART because of incremental volume updates • Instead, fold volume in half (RG) • Projection matrices are same except for reflection