Kernel-based Image Reconstruction from scattered Radon data by - PowerPoint PPT Presentation

Kernel-based Image Reconstruction from scattered Radon data by (anisotropic) positive definite functions Stefano De Marchi 1 February 5, 2016 1Joint work with A. Iske (Hamburg, D), A. Sironi (Lousanne, CH) and G. Santin (Stuttgart, D)

Main references 1 S. De Marchi, A. Iske and A. Sironi, Kernel-based Image Reconstruction from Scattered Radon Data by Positive Definite Functions , submitted 2013 (download at http://www.math.unipd.it/ ∼ demarchi/papers/Kernel based image reconstruction.pdf) 2 S. De Marchi, A. Iske and G. Santin, Kernel-based Image Reconstruction from scattered Radon data by anisotropic positive definite functions , Draft 2016 3 T. G. Feeman, The mathematics of medical imaging: a beginners guide , Springer 2010. 4 A. Iske, Reconstruction of functions from generalized Hermite-Birkhoff data . In: Approximation Theory VIII, Vol. 1: Approximation and Interpolation, C.K. Chui and L.L. Schumaker (eds.), World Scientific, Singapore, 1995, 257–264. 5 F. Natterer: The Mathematics of Computerized Tomography. Classics in Applied Mathematics, vol. 32. SIAM, Philadelphia, 2001 6 Amos Sironi, Medical image reconstruction using kernel-based methods , Master’s thesis, University of Padua, 2011, arXiv:1111.5844v1. 7 Davide Poggiali, Reconstruction of medical images from Radon data in trasmission and emission tomography , Master’s thesis, University of Padua, 2012. 8 Maria Angela Narduzzo, A kernel method for CT reconstruction: a fast implementation using circulant matrices , Master’s thesis, University of Padua, Dec. 2013. 9 Silvia Guglielmo, Stable kernel based methods for medical image reconstruction , Master’s thesis, University of Padua, Dec. 2013. 2 of 69

Part I The problem and the first approach Work with A. Iske, A. Sironi 3 of 69

Outline Image Reconstruction from CT 1 Radon transform 2 3 Alg. Rec. Tech. (ART), Kernel approach Regularization Numerical results 4 of 69

Description of CT How does it work? Non-invasive medical procedure (X-ray equipment). X-ray beam is assumed to be: - monochromatic; - zero-wide; - not subject to diffraction or refraction. Produce cross-sectional images. Transmission tomography (emissive tomography, like PET and SPECT, are not considererd here) 5 of 69

Description of CT How does it work? ℓ ( t ,θ ) −→ line along which the X-ray is moving; ( t , θ ) ∈ R × [ 0 , π ) −→ polar coordinates of line-points; f −→ attenuation coefficient of the body; I −→ intensity of the X-ray. 6 of 69

X-rays Discovered by Wihelm Conrad R¨ ontgen in 1895 Wavelength in the range [ 0 . 01 , 10 ] × 10 − 9 m Attenuation coefficient: A ( x ) ≈ ”# pho . s absorbed / 1 mm ” A : Ω → [ 0 , ∞ ) Figure: First X-ray image: Frau R¨ ontgen left hand. 7 of 69

CT machine and people Computerized Tomography (CT) modern CT Allan Mcleod Cormack Godfrey Newbold Hounsfield both got Nobel Price for Medicine and Physiology in 1979 8 of 69

Computerized Axial Tomography A. Cormack and G. Hounsfield 1970 Reconstruction from X-ray images taken from 160 or more beams at each of 180 directions Beer’s law (loss of intensity): � x 1 � I 0 � A ( x ) dx = ln I 1 Figure: First generation of CT scanner x 0 � � �� � � design. given by CT 9 of 69

Outline Image Reconstruction from CT 1 Radon transform 2 3 Alg. Rec. Tech. (ART), Kernel approach Regularization Numerical results 10 of 69

Lines in the plane A line l in the plane, perpendicular to the unit vector n θ = ( cos θ, sin θ ) and passing through the point p = ( t cos θ, t sin θ ) = t n θ , can be characterized (by the polar coordinates t ∈ R , θ ∈ [ 0 , π ) ), i.e. l = l t ,θ l t ,θ = { x := ( t cos θ − s sin θ, t sin θ + s cos θ ) = ( x 1 ( s ) , x 2 ( s )) s ∈ R } t, θ Figure: A line in the plane. 11 of 69

Radon transform definition The Radon transform of a given function f : Ω ⊂ R 2 → R is defined for each pair of real number ( t , θ ) , as line integral � � Rf ( t , θ ) = f ( x ) d x = f ( x 1 ( s ) , x 2 ( s )) ds l t ,θ R t Figure: Left: image. Right: its Radon transform ( sinogram ) 12 of 69

Radon tranform Image reconstruction A CT scan measures the X-ray projections through the object, producing a sinogram, which is effectively the Radon transform of the attenuation coefficient function f in the ( t , θ ) -plane. 13 of 69

Radon transform: another example Figure: Shepp-Logan phantom. Figure: Radon transform ( sinogram ). 14 of 69

Back projection First attempt to recover f from Rf The back projection of the function h ≡ h ( t , θ ) is the transform � π Bh ( x ) = 1 h ( x 1 cos θ + x 2 sin θ, θ ) d θ π 0 i.e. the average of h over the angular variable θ , where t = x 1 cos θ + x 2 sin θ = x T n θ . Figure: Back projection of the Radon transform. 15 of 69

Important theorems Theorem ( Central Slice Theorem) For any suitable function f defined on the plane and all real numbers r , θ F 2 f ( r cos θ, r sin θ ) = F ( Rf )( r , θ ) . (F 2 and F are the 2 -d and 1 -d Fourier transforms, resp.). Theorem ( The Filtered Back-Projection Formula) For a suitable function f defined in the plane f ( x ) = 1 2 B { F − 1 [ | r | F ( Rf )( r , θ ))] } ( x ) , x ∈ R 2 . 16 of 69

Fundamental question Fundamental question of image reconstruction. Is it possible to reconstruct a function f starting from its Radon transform Rf ? � Answer (Radon 1917). Yes, we can if we know the value of the Radon transform for all r , θ. � 17 of 69

Discrete problem Ideal case Rf ( t , θ ) known for all t , θ Infinite precision No noise Real case Rf ( t , θ ) known only on a finite set { ( t j , θ k ) } j , k Finite precision Noise in the data 18 of 69

Fourier-based approach Sampling: Rf ( t , θ ) → R D f ( jd , k π/ N ) Discrete transform: e.g. N − 1 B D h ( x ) = 1 � h ( x cos ( k π/ N ) + y sin ( k π/ N ) , k π/ N ) N k = 0 Filtering (low-pass): | r | = F φ ( r ) , with φ band-limited function Interpolation: { f k : k ∈ N } → If ( x ) , x ∈ R 19 of 69

Discrete problem Filtered Back-Projection Formula f ( x ) = 1 2 B { F − 1 [ | r | · F ( Rf ( r , θ ))] } ( x ) Filtering f ( x ) = 1 2 B { F − 1 [ F ( φ ( r )) · F ( Rf ( r , θ ))] } ( x ) = = 1 2 B { F − 1 [ F ( φ ∗ Rf ( r , θ ))] } ( x ) = 1 2 B [ φ ∗ Rf ( r , θ )]( x ) Sampling and interpolation 2 ) = 1 f ( x m 1 , x n 2 B D I [ φ ∗ R D f ( r j , θ k )]( x m 1 , x n 2 ) 20 of 69

Discrete problem: an example Figure: Reconstruction with linear interpolation and 180 x 101 = 18180 Figure: Shepp-Logan phantom. samples. 21 of 69

Outline Image Reconstruction from CT 1 Radon transform 2 3 Alg. Rec. Tech. (ART), Kernel approach Regularization Numerical results 22 of 69

Algebraic Reconstruction Techniques (ART) Differently from Fourier-based reconstruction, we consider G = span { g j , j = 1 , ..., n } of n basis functions and we solve the reconstruction problem on all Radom lines L R L ( g ) = R L ( f ) by using n � g = c j g j . j = 1 Asking interpolation, that is Rg ( t k , θ k ) = Rf ( t k , θ k ) , k = 1 , . . . , m we obtain the linear system A c = b with A k , j = Rg j ( t k , θ k ) , k = 1 , . . . , m , j = 1 , . . . , n . Large, often sparse, linear system Solution by iterative methods (Kaczmarz, MLEM, OSEM, LSCG), or SIRT techniques (see AIRtools by Hansen &Hansen 2012). 23 of 69

ART reconstruction: Example 1 Figure: 64 × 64 = 4096 reconstructed image with 4050 Figure: Bull’s eye phantom. samples by Kaczmarz. 24 of 69

ART reconstruction: Example 2 Figure: The phantom reconstructed Figure: Shepp-Logan phantom. by MLEM in 50 iterations. 25 of 69

Hermite-Birkhoff interpolation Let Λ = { λ 1 , . . . , λ n } be a set of linearly independent linear functionals and f Λ = ( λ 1 ( f ) , . . . , λ n ( f )) T ∈ R n . The solution of a general H-B reconstruction problem: H-B reconstruction problem find g such that g Λ = f Λ or λ k ( g ) = λ k ( f ) , k = 1 , . . . , n . 26 of 69

Hermite-Birkhoff interpolation Let Λ = { λ 1 , . . . , λ n } be a set of linearly independent linear functionals and f Λ = ( λ 1 ( f ) , . . . , λ n ( f )) T ∈ R n . The solution of a general H-B reconstruction problem: H-B reconstruction problem find g such that g Λ = f Λ or λ k ( g ) = λ k ( f ) , k = 1 , . . . , n . In our setting the functionals are λ k := R k f = Rf ( t k , θ k ) , k = 1 , . . . , n The interpolation conditions n � c j λ k ( g j ) = λ k ( f ) , k = 1 , . . . , n j = 1 that corresponds to the linear system A c = b as before. 26 of 69

Hermite-Birkhoff interpolation Theorem (Haar-Mairhuber-Curtis) If Ω ⊂ R d , d ≥ 2 contains an interior point, there exist no Haar spaces of continuous functions except the 1-dimensional case. The well-posedness of the interpolation problem is garanteed if we no longer fix in advance the set of basis functions. Thus, the basis g j should depend on the data: g j ( x ) = λ y j ( K ( x , y )) [= R y [ K ( x , y )]( t k , θ k )] , j = 1 , . . . , n with the kernel K such that the matrix j [ λ y A = ( λ x k ( K ( x , y ))]) j , k is not singular ∀ ( t j , θ j ) 27 of 69