18-793 Image and Video Processing
Submission instructions.
Summer 2020
• Submissions are due on Thursday 09/24 at 10.00pm ET
• Please upload scans of your solution in GradeScope (via Canvas)
Homework 3
Instructions
• Please solve all non-MATLAB problems using only paper and pen, without resorting to
a computer.
• Please show all necessary steps to get the final answer. However, there is no need to be overly elaborate. Crisp and complete answers.
• For all MATLAB problems, include all code written to generate solutions.
• Please post all questions on the discussion board on the Piazza course website.
• If you feel some information is missing, you are welcome to make reasonable assumptions and proceed. Sometimes the omissions are intentional. Needless to say, only reasonable assumptions will be accepted.
1. (Some elementary radon transforms and properties) Consider an image f(x,y) and its radon transform r(θ, α).
a) If f (x, y) is rotationally-invariant, show that r(θ, α) = r(α), i.e., the radon transform is not a function of the projection angle θ.
1 −α2 b)Letr(α)=√2πe 2 .Derivef(x,y).
(Solution:)
(a) If f(x,y) is rotationally invariant, then so is its Fourier transform. Since the FT is rotationally invariant, all radial slices are the same. Therefore, via the Fourier Slice theorem, line-projections are the same immaterial of orientation of the projection. Hence, r(θ, α) = r(α).
Mathematically,
r(α,θ) = Lθ(f)(α)
= L0(Rθ(f))(α)
= L0(f)(α) since the image is rotationally invariant
FT −ω2/2 (b) r(α) ←−−→ R(ω) = e .
2
Homework 3
From the Fourier slice theorem, we have that
R(ω) = Rθ(ω) = e−ω2/2
= F(ωcosθ,ωsinθ)
Since F(ωcosθ,ωsinθ) = e−ω2/2, setting u = ωcosθ and v = ωsinθ, we get ω =
sign(u) u2 + v2 and θ = atan(v/u) and so F(u,v)=Rθ(sign(u)√u2 +v2)=e−(u2+v2)/2
Hence, f(x,y) = 1 e−(x2+y2)/2. 2π
√
2. (Scaling and Radon Transforms) Let the Radon transform of the image i1(x, y) be r1(α, θ). Let i2(x, y) = i1(ax, ay), be a scaled image.
Find an expression for r2(α,θ), the Radon transform of i2(x,y), in terms of r1. (Solution:)
∞ −∞
∞ −∞
∞ −∞
Making the change of variables τ′ = aτ =⇒ dτ′/a = dτ,
1∞
|a|
(The absolute value compensates for the limits of integration flipping when a is negative.)
3. (Translations and Radon Transforms) Let the Radon transform of an image i1(x,y) be r1(α,θ). Leti2(x,y)=i1(x−a,y−b).
Find an expression for r2(α,θ), the Radon transform of i2(x,y), in terms of r1 and other quantities.
(Solution:)
By the definition of Radon transform,
r2(α,θ) = = =
i2(α cos θ − τ sin θ, α sin θ + τ cos θ) dτ i1(a(α cos θ − τ sin θ), a(α sin θ + τ cos θ)) dτ i1(aα cos θ − aτ sin θ, aα sin θ + aτ cos θ) dτ
r2(α,θ)= |a|
= 1 r1(aα,θ)
i1(aαcosθ−τ′sinθ,aαsinθ+τ′cosθ)dτ′
−∞
Homework 3
3
=
α−(acosθ+bsinθ)
dτ
′
α
r2(α,θ) =
= i1 Rθ τ − b
i2 Rθ τ dτ
α a
τ
= i1 Rθ τ −RθRθ b dτ τ
dτ α Ta
τ
i1 Rθ τ −(−asinθ+bcosθ) Letτ′ =τ−(−asinθ+bcosθ),anddτ′ =dτ. Thus,
τ
α−(acosθ+bsinθ) Rθ τ′
τ
i1
dτ
r2(α,θ) =
= r1(α − (a cos θ + b sin θ), θ)
4. (Implement filtered backprojection) In hw03.mat, you are given radon transform mea- surements in the variable rad, corresponding to values in alpha and theta. Implement filtered backprojection.
Notes:
(a) You are restricted to basic commands like fft, fft2, conv, conv2, meshgrid,
interp2, and other basic commands.
(b) Follow the steps laid out in the lecture: (1) 1D Fourier transform of line integral, (2)
Ramp filtering, (3) Inverse FT, and finally (4) Backprojection. (c) Attend recitation
(Solution:)
function outputImage = backproject(R, theta, recSize)
% Function to reconstruct image from Radon transform with back projection.
%
% Inputs:
% R: (Filtered) radon transform of image
% theta: Angles over which the Radon transform has been computed.
% recSize: Size of the recovered image [height, width].
%
% Outputs:
4 Homework 3
% outputImage: Recovered image with back projection.
outputImage = zeros(recSize);
% Make coordinates
[X, Y] = meshgrid(0 : recSize(1) – 1, 0 : recSize(2) – 1);
% Move origin to center of image
X = X – (recSize(1) – 1) / 2;
Y = Y – (recSize(2) – 1) / 2;
thetaRad = deg2rad(theta);
for idx = 1:length(theta)
% Smear the slice vertically
sliceProjected = ones(recSize(1), recSize(2)) .* R(:, idx)’;
% Rotate it by theta counterclockwise
Xrot = X * cos(thetaRad(idx)) – Y * sin(thetaRad(idx));
Yrot = X * sin(thetaRad(idx)) + Y * cos(thetaRad(idx));
projectionRotated = interp2(X, Y, sliceProjected, Xrot, Yrot, “linear”, 0);
% Add it to image
outputImage = outputImage + projectionRotated;
end
end
load(“hw3.mat”);
% Reconstruct size for back projection
recSize = [size(rad, 1) size(rad, 1)];
RFiltered = zeros(size(rad));
% Padding to make sure circular convolution and linear convolution give same
% results. Also powers of 2 make ffts faster.
npad = 2 ^ (ceil(log2(recSize(1)))) – recSize(1);
for idx = 1:length(theta)
% Step 1: 1D Fourier Transform of the line integral
RSlice = rad(:, idx);
RSlice = [RSlice; zeros(npad, 1)]’;
FSlice = fftshift(fft(RSlice));
Homework 3 5
% Step 2: Ramp filtering
ramp = -(length(FSlice)-1)/2 : (length(FSlice)-1)/2;
FSlice = FSlice .* abs(ramp);
% Step 3: Inverse FT, but remember to use ifftshift first
RSliceFiltered = real(ifft(ifftshift(FSlice)));
RFiltered(:, idx) = RSliceFiltered(1:size(rad, 1));
end
% Step 4: Backprojection
fbp_img = backproject(RFiltered, theta, recSize);
imshow(fbp_img, [])
The reconstructed image is shown in Figure 1.
Figure 1: Reconstructed image
Problems below wont be graded. We wont release solutions as well. We are happy to verify yours, if you post on piazza, and discuss them
5. Suppose that image f (x, y) has radon transform rf (α, θ), and the image g(x, y) has radon
transform rg (α, θ). Suppose that rg (α, θ) = d rf (α, θ). Derive an expression relating dθ
g(x,y) to f(x,y).
6
6.
7.
Homework 3
Hint: Recall that the derivative of a function a(x) at x = x0 is defined as d a(x)|x=x0 = lim a(x + h) − a(x)
dx h→0 h
Suppose that image f (x, y) has radon transform rf (α, θ), and the image g(x, y) has radon
transform rg (α, θ). Suppose that g(x, y) = d f (x, y). Derive an expression relating rg in dx
terms of rf .
You are given the radon transform of an image. Derive an analytical expression for the image. It is ok if the values in your expression are approximate.
-150
-100
-50
0
50
100
150
1.5
1
0.5
0
0 20 40 60 80 100 120 140 160 theta in [degrees]
alpha