Magic Inverse Problems Course - The Radon Transform Bill Lionheart - PowerPoint PPT Presentation

Magic Inverse Problems Course - The Radon Transform Bill Lionheart March 25, 2013

The line L Θ , s with normal vector Θ = (cos θ, sin θ ) and distance s to the origin is the set L Θ , s = { x ∈ R 2 | x · Θ = s } . Θ takes its values on a circle S 1 = { (cos θ, sin θ ) | 0 ≤ θ < 2 π } so the space of lines in R 2 is S 1 × (0 , ∞ ) (or maybe S 1 × R but that counts every line with s � = 0 twice). Topologically this is a cylinder. We define the Radon transform of a function f : R 2 → R to be a function R [ f ](Θ , s ) on the set of lines � � R [ f ](Θ , s ) = fdl = fdl . L Θ , s x · Θ= s We can parameterize a line using the coordinate t as distance along the line. Define Θ ⊥ = ( − sin θ, cos θ ), unit vector perpendicular to Θ. Then any x point on L Θ , s has the form x = s Θ + t Θ ⊥ , in this form � � ∞ � s Θ + t Θ ⊥ � fdl = f dt . x · Θ= s t = −∞

Fourier Slice Theorem We can use the Fourier Transform to better understand the Radon transform. We need the F-T. of a function of 2 variables—or in general n -variables � 1 R n f ( x ) e − i x · w d x . ˆ f ( w ) = (2 π ) n / 2 Some authors miss out the constant, but with this definition, the inverse F-T. is � 1 R n f ( w ) e i x · w d x ˇ f ( x ) = (2 π ) n / 2 and of course ˇ f = ˆ ˆ ˇ f = f . Also the scaling, translation, differentiation and Perseval’s properties hold from the case n = 1. Let us define, for a fixed Θ R Θ [ f ]( s ) = R [ f ](Θ , s ) for any s ∈ R . This is all the lines in a particular direction we have.

Theorem (Fourier Slice Theorem) √ � 2 π ˆ R Θ [ f ]( σ ) = f ( σ Θ) . Here σ ∈ R is the F-T.(frequency) variable for s. Proof. � ∞ � 1 f ( x ) e − is σ dl ds � R Θ [ f ]( σ ) = √ 2 π −∞ Θ · x = s � ∞ � ∞ 1 f ( s Θ + t Θ ⊥ ) e − is σ dt ds . = √ 2 π s = −∞ t = −∞ Now � x � � � � s � cos θ − sin θ x = = sin θ cos θ y t which is simply a rotation of the coordinate system ( x , y ) to ( s , t ). � � � � ∂ ( x , y ) � � � � � = 1 . � ∂ ( s , t ) � ∞ � ∞ 1 � f ( x ) e − i σ s dt ds R Θ [ f ]( σ ) = √ 2 π −∞ −∞ � ∞ � ∞ 1 f ( x , y ) e − i σ Θ · x dx dy = √ 2 π −∞ −∞ √ 2 π ˆ = f ( σ Θ) .

So by taking the Radon transform data, Fourier transformed in s , gives us the Fourier transform of the image along each direction. At least where the F.Ts exist in the classical sense this shows that we have uniqueness of solution for the inverse Radon transform. Corollary � R 2 � For f ∈ S , R [ f ] uniquely determines f (that is R is injective). Proof. From FS Theorem R [ f ] determines ˆ f on R 2 , hence f by the Fourier inversion formula. Note, althoughˆ: L 2 � R 2 � → L 2 � R 2 � is invertible, R isn’t necessarily defined on L 2 � R 2 � , as as a line is a set of measure zero. You need more smoothness for line integrals to be defined. It is possible to prove for example that R is defined and continuous, with a continuous inverse in the Sobolev spaces � R 2 � and H s +1 / 2 (lines)(see [1] p. 42). (Not quite honest, I should say defined on a bounded subset of R 2 ). H s 0

Backprojection The Fourier Slice Theorem is not as useful as you might first think. The problem is that if we discretize Θ and σ in regular steps, this produces a non-uniform image resolution. To derive other reconstruction algorithms we must consider a kind of adjoint. Recall that for K : H 1 → H 2 , � K [ f ] g � = � f , K ∗ [ g ] � defines the adjoint. For the special case of R n and R m , K could be a matrix and K ∗ its transpose. The ‘formal adjoint’ of our operator K defined on functions is the adjoint with respect to the L 2 norm, even if K isn’t actually defined on L 2 . That is for all f and g for which the integrals make sense � � f K ∗ [ g ] , K [ f ] g = where K ∗ is the formal adjoint, for example let K [ f ] = f ′ on functions on R vanishing outside some interval, then � ∞ � ∞ � ∞ � ∞ ′ f ′ g = [ fg ] ∞ K [ f ] g = −∞ − fg = − fK [ g ] . −∞ −∞ −∞ ∞ So K ∗ = − K , but it isn’t K ∗ as it’s not defined on L 2 ( R ). The formal adjoint of R , is an operator R ∗ which takes Sinogram data, functions on S 1 × R and turns it into ‘images’, functions on R 2 . � � � R 2 R ∗ [ g ]( x ) f ( x ) dx = S 1 R [ f ](Θ , s ) g (Θ , s ) . R So clearly R ∗ [ g ]( x ) is the integral of g over all lines through x . It’s called ‘backprojection’ as it throws the data back into the image space. It doesn’t reconstruct however.

In more specific terms g (Θ , s ) backprojected � 2 π R ∗ [ g ]( x ) = g (Θ , Θ · x ) d θ, θ =0 where Θ = (cos θ, sin θ ) as usual, and to make the integral over the (Θ , s ) pairs with Θ · x = s , we have substituted for s . Visually, the effect of R ∗ R on a small blob is to smear the value of a line integral along a particular line through the blob, all over the line. This turns a blob into a star, at least when infinitely many angles are used. If one attempted to use ‘backprojection’ as a reconstruction algorithm one would get a ‘starry’ version of the image. This would be the first iteration of Landweber’s method for solving R [ f ] = g , starting with f 0 = 0. Landweber like methods have been used in X-Ray CT where they are called SIRT—simultaneous iterative reconstruction technique, as all the data is ‘backprojected’ at once. By construct ART—algebraic reconstruction technique ‘backprojects’ one block. e.g., R Θ at a time. This is Karczmarz.

Filtered backprojection This is done in more detail in [1], p. 18–23. We simplify the argument by taking n = 2. We define the ‘ramp filter’ operator on a function f : R → R by H [ f ]( w ) = | w | ˆ � f ( w ) (It is not standard notation, Natterer calls it I − 1 ). f , so H ◦ H = − ∂ 2 /∂λ 2 the second Note that it is a bit like a differential operator, indeed � H ◦ H [ f ]( w ) = | w | 2 ˆ � − ∂ 2 /∂λ 2 . derivative So we could legitimately say H = Now by definition of the I.F-T. � H [ f ]( x ) = (2 π ) − 1 / 2 R 2 e i x · w | w | ˆ f ( w ) d w . (notice here f is a function of one variable) Using polar coordinates in the w plane w = σ Θ and the Fourier inversion formula. � 1 R 2 e i x · w ˆ f ( x ) = f ( w ) d w 2 π � 2 π � ∞ 1 e i x · σ Θ ˆ = f ( σ Θ) σ d σ d θ. 2 π θ =0 σ =0

The Fourier Slice Theorem says √ 2 π ˆ f ( σ Θ) = R Θ [ f ]( σ ) so � 2 π � ∞ 1 e i x · σ Θ R Θ [ f ]( σ ) σ d σ d θ f ( x ) = (2 π ) 3 / 2 θ =0 σ =0 � 2 π � ∞ 1 e i x · σ Θ σ ˆ = g (Θ , σ ) d σ d θ, (2 π ) 3 / 2 θ =0 σ =0 where g = R [ f ] and ˆ g is the F-T. w.r.t the second variable s . � 2 π � ∞ 1 1 e i x · σ Θ | σ | ˆ f ( x ) = g (Θ , σ ) d σ d θ (2 π ) 3 / 2 2 θ =0 −∞ (remembering the symmetry of g ) � 2 π 1 1 = H [ g ](Θ , Θ · x ) d θ 2 2 π θ =0 1 R ∗ H [ g ] . = 4 π So in principle we can reconstruct Radon transform data by ‘ramp filtering’ the data then ‘backprojecting’.

More general reconstruction In the more general case (still n = 2) Natterer defines I α f ( w ) = | w | − α ˆ � f ( w ) . So H = I − 1 which is defined on data (on the cylinder) using the FT wrt the second variable, and on ‘images’ R 2 → R 2 using the 2D.FT. He then shows 1 (2 π ) − 1 I − α R ∗ I α − 1 g , f = 2 where g = Rf , which is a reconstruction using a filter on the data and on the backprojected image. For α > 0, the image is smoothed (high frequencies suppressed) and high frequencies in the data are sharpened. This type of technique can result in a practical reconstruction algorithm. (Note that on S , I α I − α is the identity.) 1 (2 π ) − 1 I − 1 R ∗ R [ f ] f = α = 1 . 2 So R ∗ R [ f ] = 4 π I 1 [ f ], from which we can see that the Tikhonov formula � R ∗ R + µ 2 � − 1 R ∗ g f µ = can be implemented as a filter applied to backprojected data, or conversely f µ = R ∗ � RR ∗ + µ 2 � − 1 g backprojecting filtered data.

SVD of Radon Transform See [2] p. 240–245, and [1]p. 95–101. In some ways Natterer’s treatment is clearer—but he does it in n -dimensions. We will consider our ‘image’ f to have support contained in a unit disc centered on the origin D , and f ∈ L 2 ( D ). The codomain of the Radon transform, functions on the cylinder, will need a certain weighted norm. � 2 π � 1 | g ( s , Θ) | 2 � g � 2 y = � 1 − s 2 ds d θ. 0 − 1 We will call this weighted L 2 space Y and have X = L 2 ( D ), R : X → Y . (Natterer shows that R and R ∗ are both bounded on these spaces.) Let’s check this. � Define w ( s ) = 1 − s 2 � w ( s ) � s Θ + t Θ ⊥ � R [ f ]( s , Θ) = f dt | s | < 1 − w ( s ) � w ( s ) � � s Θ + t Θ ⊥ �� 2 dt . � � | R [ f ]( s , Θ) | 2 ≤ 2 w ( s ) � f � − w ( s ) So � 1 � 1 � w � � s θ + t Θ ⊥ �� 2 dt ds � � ( w ( s )) − 1 | R [ f ]( s , Θ) | 2 ≤ 2 � f � − 1 − 1 − w � | f ( x ) | 2 dx = 2 D � 2 π � 1 | g ( s , Θ) | 2 � g � 2 y = ds d θ. w ( s ) 0 − 1 So we have explicitly √ � R [ f ] � y ≤ 4 π � f � x .