Medical Image Analysis
Medical Image Analysis CS 778
Popular in Course
Popular in ComputerScienence
This 82 page Class Notes was uploaded by Abe Jones on Saturday September 12, 2015. The Class Notes belongs to CS 778 at West Virginia University taught by Timothy McGraw in Fall. Since its upload, it has received 43 views. For similar materials see /class/202761/cs-778-west-virginia-university in ComputerScienence at West Virginia University.
Reviews for Medical Image Analysis
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 09/12/15
Medical Image Analysis CS 778 578 Computer Science and Electrical Engineering Dept West Virginia University January 28 2009 39 Outline 0 Anisotropic diffusion 9 Energy Minimization 39 Outline Anisotropic diffusion o PeronaMalik 0 Edge enhancing 0 Coherence enhancing 0 Energy Minimization 2059 3 m3 PezonaeMahk Perona Malik Implementation 8M diV GEM CW1 u csuy 7 CE 0 CW 0 7 d1v 0 EN Vu 0 CS V14 The implementation is anisotropic diffusion with a diagonal diffusion tensor JanuarylELZOGQ 2 23 Edge39exilmnclug The diffusion tensor D From the statistical mechanics of diffusion it is known that D is the covariance matrix of the molecular displacement probability density function The probability that a molecule at position x0 at time t will diffuse to x0 r at time t 739 p7 r N07 27D So D must be 0 Symmetric D DT a Nonnegative de nite xTDx 2 0 for all Xy 0 January2amp 2009 5 23 Edge39enlmnclug The diffusion tensor D Physically these constraints mean 0 Symmetric conservation of mass 0 Nonnegative de nite no backward diffusion VMTDVM gt 0 We can understand the geometry of anisotropic diffusion by looking at the eigenvalue decomposition of D D XAX l JanuarylELZOO9 5 23 l Edge39enlmnclug The diffusion tensor D D XAX l o X el 82 e are the eigenvectors ofD 7 A10 Ailo 2l 0 A1 2 are the eigenvalues of D 0 Since D is real and symmetric the eigenvalues are real 0 Since D is symmetric the eigenvectors are mutually orthogonal 0 Since D is nonnegativede nite the eigenvalues 2 0 January28gl 9 7 x23 l gamma Eigenvalue Decomposition of D I SinceX is orthogonal X 1 XT a X is a rotation matrix 0 A is a nonuniform scaling matrix D XAXT l Edge39enlmnclug Eigenvalue Decomposition of D 0 Since X is orthogonal X 1 XT o X is a rotation matrix 0 A is a nonuniform scaling matrix D XAxT o Rotate ux by multiplying with rotation matrix 0 Perform nonuniform scaling 0 Apply inverse rotation JanuarylE OOQ gjquotzs Eageren xmaug Construct D Construct X using regularized smoothed edge information Let MU K0 u 0V1 HVMU OVZLVI HV1HHV2H1 0XV1 V2 7 1 0 Ho A2 0 1 g Wad the diffusivity from PeronaMalik 0 2 1 Januanyizmoogg 110923 Eagere imoiug Construct D When g z 1 away from an edge then D x I andj z 7Vu isotropic diffusion Nearanedgeg a 0 Eagee lmuclug Construct D DWV1V2lt3 H MW DVMV1 V2V29VM DVM V2V2 V14 DVu is parallel to V2 so it is parallel to the edge Edge enhancing Tensor eld for input image Note that tensors are isotropic Within homogeneous regions Tensors are anisotropic and direct the ux along edges The resulting diffusion process is edge enhancing es 727 mags margarita ft f g i January 28 2009 13 23 39 Coherence enhancing diffusion Structure matrix JPVMU KP VMUVug 172 M M J V14 1 Ly plt 0 Let V1 and V2 be the eigenvectors of J and p1 and M be the corresponding eigenvalues Jamammo 39 Structure matrix a Direction V1 has the greatest intensity variation 0 Direction V2 has the least intensity variation 0 M m give the degree of intensity variation contrast in the eigendirection o Constantvalued image a 11 p2 0 0 Step edge a p1gt 07 z 0 o Corneralul gt 07M gt 0 January 283009 15 23 l In computer Vision D has the same eigenvectors as the structure matrix but different eigenvalues A1 A2 A1204 oz ifm M2 2 oz l oz exam otherwise When m LL2 are very different smooth more in the coherence direction CS 778 5 78 West Virginia Universityquot Medical Image Analysis 39 Outline 0 Anisotropic diffusion 9 Energy Minimization o Variational calculus 0 Imposing data constraints 39 Variational Calculus There is another way to arrive at the diffusion equation We wish to nd u such that E is minimized where E is a functional of the form Eltugt Mmmunuymx A functional depends not just on variable Xy but also on functions and their derivatives 14 Ml My We will work with functionals by treating 147141714 as variables January 28mm lg 23 39 Variational Calculus When E is a functional of the fonn Eltugt fltx7y7u7unuygtdx Q then when E is minimized the following condition holds 8 8 WE fu guxiaiyfuy0 fu 7 d1V f fay The minimum can be found using the descent method 8M 7f diV g January 283009 20 23 39 Variational Calculus Membrane spline energy represents the stretching arclength of a thin sheet Minimizing this energy results in smooth surfaces The membrane spline energy has the form 19mmAHWHMAM M x January 283009 3921 23 39 Variational Calculus Minimizing the membrane spline energy m HVMHde 14 54x 0 n has EulerLagrange condition for minimization 8 8 WE faiawgy yd 8 8 0 7 gal197 gal1y 0 Leading to the descent equation 8M diVVu January 283009 22723 39 Variational Calculus Later we will look at a more general form of the variational calculus which will allow us to minimize other functionals such as the thinplate spline energy This functional minimizes bending energy curvature ETPSM u 214iy a dx 9 January 283009 23 23 l Data Constraints Let M0 be the original input image We wish to penalize solutions which are far from the initial condition Find min Eu m u 145 u 7 m2 dx 9 2 results in the descent equation a reactiondiffusion equation 8M diVVu M140 7 u The evolution reaches a nontrivial steadystate January 28mm 24 23 Sparse block matrices in Matlab Constructing sparse block matrices Sparse matrices Block matrices Solving large sparse linear systems LU factorization Conjugate gradient Sparse matrices sparsemn All zero sparse mxn matrix sparseA Converts full matrix A to sparse Speyemn Sparse matrix with ones on the main diagonal spallocmnnz Allocates storage for an mxn matrix with nz nonzero entries Since reallocation is expensive it is a good idea to allocate storage for a matrix before building it Sparse matrices spdiagsB d m n Form a sparse mxn matrix whose diagonals d are the columns of B lnd O is the main diagonal Positive values are super diagonals Negative values are subdiagonals Example second central difference matrix e ones41 A spdiagse 2e e 1 0 1 4 4 2 1 0 A 1 2 1 0 0 1 2 1 0 1 2 Sparse matrices 39 SPYA Visualize the sparsity structure of the matrix Example 1D second central difference matrix n 32 e onesn1 A spdiagse 2e e 1 O 1 n n SpyA Block Matrices Sometimes it is useful to specify a matrix blockbyblock M blkdiagab a11 a1n O O M an an O O O O b11 b1 blkdiag example e ones33 A blkdiageee spyA Matrix concatenation horzcata1 a2 a3 Concatenate matrices horizontally vertcata1 a2 a3 Concatenate matrices vertically Matrix concatenation e ones33 39 SPYA z zeros33 eye33 A horzcate Z l 1 2 3 A 5 E 7 E a 1 quot212 B vertcate z I Kronecker tensor product K kronXY ifX is mxn and Y is pxq then K is mp x nq Kronecker tensor product example spykronXY X ones33 Y eye33 39 SPYkrOHY X Using kron to create a 2D Laplacian matrix Boundary conditionszzeros outside image domain First create 1D second central difference matrix for x direction n1 size1 e1 onesn11 1 speyen1 n1 D1xx spdiagse12e1e11 01 n1 n1 SpyD1XX 1 2 3 A 5 quot213 Then create the 2D second central difference matrix I2 speyen2 n2 D2xx kron2 D1xx spyD2xx u u u 5 u u u u u m u u u u u 15 u u u cu u 2n n u u u u 25 5 m 15 2m 25 quot255 Using krcn to create a 2D Laplacian matrix Create 1D second central difference matrix for ydirection n2 size1 e2 onesn21 I2 speyen2 n2 D1yy spdiagse2 2e2 e2 1 0 1 n2 n2 Then create the 2D second central difference matrix D2yy kronD1yy l1 spyD2yy 2D Laplacian Matrix Compute 2D Laplacian matrix L D2xxD2yy u u u 5 5 u 5 u u u 1U 1n n 1n n u u u 15 15 u 15 u u u u 2U 2U cc 2n n u u u 25 25 u 25 5 1U 15 2D 25 quot255 In 3D D3xx kron3 kron2 D1xx D3yy kron3 kronD1yy 1 D322 kronkronD1zz I2 I1 L D3xxD3yyD322 Imposing other boundary conditions Periodic boundary conditions D1xx D1XX spdiagse1e1 n11 n1 1 n1 n1 D2XX kronl2 D1XX D1yy D1yy spdiagse2 e2 n21 n2 1 n2 n2 D2yy kronD1yy I1 Solving linear systems Solve forx in Axb Inversion xA4b Not practical for large or illconditioned matrices Other direct methods LU factorization Iterative methods Conjugate gradient CG methods LU factorization This may be what happens when you type 39x Ab39 in Matlab Check mldivide help for details LU decomposition is a form of Gaussian elimination Permits the linear system to be solved by back substitution If the matrixA does not change in every iteration you can factorize the matrix once then only perform the back substitution each iteration LU factorization A LU L is lower triangular all superdiagonals are 0 U is upper triangular all subdiagonals are 0 11 0 0 0 11 12 13 um 21 22 0 0 0 22 23 2n L 31 32 33 0 U 0 0 33 3n n1 n2 n3 rm 0 O 0 ml For symmetricA you can find the Cholesky decomposition ALU Solving by LU factorization ReplaceAwith L times U Axb LUxb Solve in 2 steps Let yzUx Solve Lyb Then solve Ux y Solving triangular linear systems Easy just back substitution Proceed rowbyrow Solve for one unknown per row 11 0 0 0 3 1 b1 21 22 0 0 3 2 b2 Lzl b 31 32 33 0 3 3 3 2 Solving by LU decomposition in Matlab To solve Ax b Decompose LU MN Backsubstitute x ULb Sparse LU lfA is sparse then L and U are usually sparse also For the 2D Laplacian matrix 5 5 m m L U m m m m K K 5 l 15 2D 25 5 lEI 15 2D quot2 EU quot2 2W Conjugate gradient Iterative method Only requires matrixvector multiplications vector vector operations Can solve symmetric positivedefinite systems See JR Shewchuk An introduction to the conjugate gradient method without the agonizing pain for more details Conjugate gradient Green lines iterations of gradient descent Subsequent search directions v are perpendicular 39 ViTVi1 0 Red lines iterations of conjugate gradient method In CG methods the search directions are conjugate viTAvi1 O Conjugate gradient variants in Matlab Preconditioned CG symmetric A x pcgAbtomaxitM Bico njugate gradients square A not req39d to be symmetric X bicgAbtomaxitM I CG squared a variant of bicg x cgsAbtomaxitM Biconjugate gradients stabilized method anothervariant of bicg x bicgstabAbtolmaxitM See Matlab help for details on differences in computational cost and convergence speed Preconditioning The matrix M specified in the Matlab functions is a preconditioner M1AxM1b lfA is illconditioned choose M such that M39lA is well conditioned M must be symmetric and positive definite for PCG Ideally M391 A391 Medical Image Analysis CS 778 578 Computer Science and Electrical Engineering Dept West Virginia University March 9 2009 39 Outline 0 Image Registration 0 Classi cation of Registration Methods 9 Overview of imaging modalities 39 Outline 6 Image Registration 39 Problem De nition Image registration is the process of determining a coordinate transformation between two images that are misaligned me distI1 X 12TX o T is a coordinate transformation 0 1 and 2 are 2 images to be aligned o distIlIz is a metric which determines how well the images match 0 dist11 2 can be based on image intensities or extracted features I Maich Q 2009 4 M7 Image Registration Examples Panorama gtgfiing CS 778 578 West Virginia University March 92009 5 47 Wm Examples Multimodal rigid registration CT MR registration March 9 2009 6 M7 Examples March 9 2009 747 I Some Applications 0 Determine a correspondence between function PET fMRI and anatomy MRI CT 0 Longitudinal studies of individuals over time 0 Comparison of different subjects MafchBQZOGQ 8 m 39 Outline 9 Classi cation of Registration Methods 39 Criteria From A Survey of Medical Image Registration by J Maintz and M Viergever 0 Dimensionality of images 9 Nature of registration basis markers landmarks intensity 9 Nature of transformation translation rotation elastic 9 Domain of transformation global local 3 Interaction interactive automatic 9 Optimization procedure computation search 0 Modalities involved CT MRI 9 Subject Same patient 0 Object anatomical region Mazph 92009 1047 l 1 Dimensionality Medical images may be 2D 3D or even 4D 0 2D projections Xray or slices of volumes 0 3D volumes CT MR 0 4D time series volume images cardiac motion We may wish to register 0 a Spatial Dimensions 0 2D 2D 9 2D 3D 9 3D 3D 0 b Time Series 0 2D 2D 9 2D 3D 9 3D 3D 0 3D 3D is most common 0 Increasing dimension leads to more complex transformations and metrics Mmm zzoog 1147 l 2 Nature of registration basis 0 a Extrinsic foreign objects in the image volume 0 b Intrinsic based on the subject being imaged o c Nonimage based calibrated coordinate systems multimodal scanners PETCT hybrid ilk3947 Classi cation of Registration Methods 2 Nature of registration basis 0 a Extrinsic 0 Invasive A Stereotactic frame B Fiducials screw markers 9 Non invasive A Foam mold frame dental adapter B Fiducials skin markers v Stereotactic frame left and resulting image center Noninvasive skin markers right CS 778 578 West V aUniversity 39 u 39 7 739 March 9 2009 13 47 l 2 Nature of registration basis 0 b Intrinsic 0 Landmark based at A Anatomical B Geometrical 0 Segmentation based at A Rigid models points curves surfaces B Deformable models 9 Voxel property based A Reduction to features centroid principal axes at B Using full image content ation of Registratixm Methods 1 Landmark based Anatomical Landmarks Sparse set of points based on known anatomy Anterior commissure Posterior commissure Geometrical Landmarks Corners curvature maxima CS 778578 1W March 92009 15 Classi cation of Registration Methods 2 Segmentation based Segment the same structure from both images then align these structures Rigid Models Matching edges and surfaces Deformable Models Deformable atlases FEOX39ITAL Accuracy of the registration is limited by the accuracy of the segmentation CS 778 578 West Virginia University March 92009 16 47 138me ation of Registration Methods 3 Voxel property based Reduction to features an inn an 2m 25m 3m an inn an 2m 25m 3m an inn an 2m 25m 3m Centroid principal axes of the images Using full image content Image intensity March 9 2009 1747 ion M ethods Nature and domain of transformation Original Global Local 0 3 Nature of transformation gt a Rigid gt b Af ne gt c Projective F d Curved o 4 Domain of transformation Projective gt a Local rarely used Curved gt b Global CS 77 8 57 8 Vest Tniversityh 39 March 9 2009 18 l 47 39 Nature and domain of transformation Represented by matrices o a Rigid Rotation translation re ection o b Af ne shear o c Projection 2D3D registration Classi cation of Registration Methods Nature and domain of transformation 0 d Curved Elastic Polynomial transformation b splines CS 778 578 West Virginia University Medical Image Analysis March 9 2009 47 l 5 Interaction o a Interactive 0 Initialization supplied 9 No initialization supplied 0 b Semiautomatic 0 User initializing 9 User steering correcting 9 Both 0 0 Automatic l 6 Optimization procedure a a Parameters computed explicitly a b Parameters searched for iterative techniques Explicit computation is generally possible only when the number of parameters is small as when registering small point sets Search Iterative techniques like gradient descent may be used to optimize the energy function I 92 009 22747 l 7 Modalities involved a a Monomodal gt Time series gt Pre and post contrast enhancement 0 b Multimodal gt Structural Structural gt Structural Functional o c Modality to model Mathematical models of anatomy or pathology 0 d Patient to modality interoperative or radiotherapy applications l 8 Subject 0 a Intrasubject same patient 0 b Intersubject different patients 0 0 Atlas statistical database Classi cation of Registration Methods Atlas 0 Statistical atlases are models of expected or normal anatomy 0 May be constructed by sampling from large population of images 0 These samples are registered to a common coordinate system and averaged CS 778 578 West Virginia University March 9 2009 25 47 l 9 Object Most frequently 0 a Head 0 b Thorax o c Abdomen o d Pelvis o e Limbs o f Spine and vertebrae 39 Outline 9 Overview of imaging modalities Xray 0 Physical basis Xrays are attenuated absorbed or scattered by hard tissue and may pass through soft tissue 0 Hardware Xray tube and detector a Image formation Photographic plate or electronic detector scintillator 0 Applications Skeletal structure mammography angiography with contrast agent 92 009 28 m Xray M djcal IthagttA 39 March Q 2009 29 47