Description
1.2 Library imports and writing tests
The code we provide includes a superset of all the imports you are allowed to use when implementing your
solution.
One exception to this policy is that, in the __main__ block, you are allowed to import any additional
routines you would like in order to test your solution code. Do not use import in any of your solution
routines nor at the top of the Python file; only optionally in __main__.
2 Discrete Fourier Transform
Given an input continuous function sampled at equally-spaced discrete points on an integer lattice
in the primal domain (e.g., time, space, etc.), we denote this discretized signal as
. Note that, in general, can be a complex-valued function (and, so too the sampled values
); however, our examples will only consider real-valued input functions .
The discrete Fourier transform (DFT) of the sampled signal is then given by:
which can be equivalently expressed as a matrix-vector operation (see lecture slides) as
You must not use any additional imports for your solution, other than the ones provided by
us.
Do not edit the imports at the top of the .py file. Doing so will result in an assignment grade of
0%.
☒
It is very difficult to successfully complete this assignment without writing your own tests in
__main__. We strongly recommend against relying exclusively on the example tests we provide.
Think about what kind of input/output behaviours you can write and test against.
ⓘ
𝑓 𝑁
𝑥 ∈ {0,…, 𝑁 − 1}
𝑓(𝑥) ≡ 𝑓[𝑥] 𝑓
𝑓[⋅] 𝑓 : ℝ → ℝ
𝑓[𝜔]
̂
[𝜔] = 𝑓[𝑘] exp( )
𝑓 ,
̂ ∑
𝑘=0
𝑁−1 −2𝜋ı𝑘𝜔
𝑁
(1)
ˆ𝐟 = 𝐌𝑁 𝐟 ,
where is a vector of all the sampled function values, encodes the
linear operator form of the DFT for a length- input, and is a vector of
the sampled DFT of . The (discrete) signal can be interpreted as a discrete sampling of the
continuous Fourier transform of an -periodic summation of — i.e., repeated
periodically.
Note that, regardless of whether the input signal is real- or complex-valued, the DFT always has
complex-valued elements: encodes the frequency amplitudes and their phases.
The inverse discrete Fourier transform (iDFT) allows us to transform back from this (discrete, complexvalued) sampled frequency spectrum to the original (possibly only real-valued) sampled signal in the primal
domain, as
Note the similarities in the mathematical definitions of the DFT and iDFT in Equations and , comprising of
a sign change in the complex exponential and different normalization constants. You will leverage these
similarities to avoid code duplication: your DFT and iDFT implementations will rely on an implementation of a
generalized (i.e., bidirectional) DFT.
Specifically, you will implement this generalized DFT formula for a sampled input ,
where we obtain the DFT by setting and (and and ), and the iDFT with
and (and and ). Your implementations will act on, and return, an entire
vector of sampled primal- (i.e., ) or frequency-domain (i.e., ) values.
As discussed above, the generalized DFT equation can be used to realize both a forward and inverse DFT.
𝐟 = (𝑓[0], 𝑓[1],…, 𝑓[𝑁 − 1]) 𝐌𝑁
𝑁 ˆ𝐟 = (𝑓[0], [1],…, [𝑁 − 1])
̂ 𝑓̂ 𝑓̂
𝑓 𝑓[𝑥]
̂
{𝑓𝑁 (𝑥)} 𝑁 𝑓(𝑥) 𝑓(𝑥)
𝑓[𝑥]
̂
Re(𝑓)
̂ Im(𝑓)
̂
𝑓[𝑥] = [𝑘] exp( )
.
1
𝑁 ∑
𝑘=0
𝑁−1
𝑓̂
2𝜋ı𝑘𝑥
𝑁
(2)
1 2
Υ[⋅]
Ψ[𝑎;s] = Υ[𝑏] exp(
s ) ∑ ,
𝑏=0
𝑁−1
2𝜋ı𝑏𝑎
𝑁
Υ = 𝑓 𝑠 = −1 𝑎 ≡ 𝜔 𝑏 ≡ 𝑥
Υ =
1
𝑁
𝑓̂ 𝑠 = 1 𝑎 ≡ 𝑥 𝑏 ≡ 𝜔
𝐟 ˆ𝐟
Deliverable 1 [17.5%]
Implement the Generalized DFT Routine: Complete the implementation of the generalized DFT
function DFT as specified by and with the default parameter setting of .
☑
Ψ s = −1
Υ
Deliverable 2 [2.5%]
Implement the inverse DFT Routine: Complete the implementation of the inverse DFT function
iDFT that relies on calling your generalized DFT routine, above.
☑
3 2D Discrete Fourier Transform
Let’s now consider a generalization to 2D, where we sample a continuous function on the
(square) 2D integer lattice leading to a discretized 2D signal at .
We can similarly define the sampled 2D DFT of as
which admits a similarly compact matrix-vector expression as
where we can arrange the sampled primal function values in a matrix with elements , with
encoding the linear operator form of the 2D DFT for a square length- sampled input, and
is the matrix of the sampled 2D DFT coefficients with .
Again, similarly to the 1D setting, DFT elements are always complex-valued with encoding
frequency amplitudes and their phases.
The 2D inverse DFT is now
and a similar extension of the generalized DFT routine applies in the 2D setting.
3.1 Leveraging Factorization
After little manipulation, one can observe that the 2D DFT and iDFT in Equations and can be factorized
into several applications of their 1D counterparts.
A careful numpy implementation should not rely on any for loops or related Python operations (e.g.,
loop comprehensions). You may also wish to leverage the fact that exp (ı𝑥) = cos 𝑥 + ısin 𝑥.
𝑓 : ℝ → ℝ
2
𝑓[𝑥, 𝑦] (𝑥, 𝑦) ∈ {0,…, 𝑁 − 1}
2
𝑓[ , ]
̂𝜔𝑥 𝜔𝑦 𝑓[𝑥, 𝑦]
[ , ] = 𝑓[ , ] exp ( )
𝑓 ,
̂𝜔𝑥 𝜔𝑦
𝑘
∑𝑥=0
𝑁−1
𝑘
∑𝑦=0
𝑁−1
𝑘𝑥 𝑘𝑦
−2𝜋ı(𝑘𝑥𝜔𝑥 + 𝑘𝑦𝜔𝑦)
𝑁2
(3)
ˆ𝐅 = 𝐌𝑁 𝐅 , ,𝑁
𝐅 (𝐅)𝑖,𝑗 = 𝑓[𝑖,𝑗]
𝐌𝑁,𝑁 𝑁 × 𝑁 ˆ𝐅
(ˆ𝐅 )𝑖,𝑗 = 𝑓[𝑖,𝑗]
̂
𝑓[⋅, ⋅]
̂ Re(𝑓)
̂
Im(𝑓)
̂
𝑓[𝑥, 𝑦] = [ , ] exp ( )
1
𝑁2
𝑘
∑𝑥=0
𝑁−1
𝑘
∑𝑦=0
𝑁−1
𝑓̂𝑘𝑥 𝑘𝑦
2𝜋ı(𝑘 𝑥 + 𝑦) 𝑥 𝑘𝑦
𝑁2
(4)
3 4
By applying a 1D DFT independently to each of the rows of , before applying the 1D DFT independently to
the columns of the (partially-)transformed , we arrive at the 2D DFT transform of . In matrix form, this
amounts to
You should leverage the factorization to accelerate your implementation.
4 Fourier for Image Convolution & Deconvolution
In Assignments 1 and 3, we explored convolution and deconvolution using linear operators that act on
sampled signals in the primal (i.e., pixel, spatial) domain.
We constructed the convolution operators as a function of the underlying convolution kernel, and we
implemented a wrap-based boundary condition in instances during convolution where kernels only partially
overlapped the input domain.
Manually treating the indexing, which additionally required “flattening” 2D images into 1D vectors, required
care.
Frequency-space representations are particularly well-suited to certain applications, among them
convolution. We will use the DFT to perform image convolution and (regularized) deconvolution, leveraging
these advantages:
1. convolution in the primal domain amounts to multiplication in the frequency domain, and
2. in the frequency domain, we can perform computation directly with native 2D signals (e.g., the image
and kernel), instead of having to “flatten” them for the sake of formulating the forward (and backward)
problems as systems of linear equations.
𝐅
𝐅 𝐅
ˆ𝐅 = 𝐌𝑁 𝐅 = ,
,𝑁 𝐌𝑁 (𝐌𝑁 𝐅 )
T T
(5)
Deliverable 3 [22.5%]
Implement a Generalized 2D DFT Routine: Complete the implementation of the generalized 2D
DFT function DFT2D, which extends the 1D generalized routine to 2D.
☑
Keep in mind that you don’t have to explicitly form these matrices at any point in order to take
advantage of the aforementioned factorization property.
Deliverable 4 [2.5%]
Implement the inverse 2D DFT Routine: Complete the implementation of the inverse 2D DFT
function iDFT2D that relies on calling your generalized DFT2D routine.
☑
4.1 2D Fourier Convolution
The (forward) convolution of a 2D function with a 2D kernel yields a 2D output function
. In the image convolution settings we have explored, is a
discretely-sampled input image, is a (typically square, odd edge-length) discrete kernel,
and is the blurred output image.
Given a matrix of the 2D sampled image pixel values, i.e., with , we obtain its DFT as per
Equation . For the kernel, we can also form a matrix of its sampled values before computing its DFT ;
here, however, we need take into account two additional points:
1. the kernel matrix needs to be 0-padded to match the size of the image matrix — this will be
needed to allow for an element-wise product of their DFTs, later on; and,
2. the “location”/placement of the non-zero kernel elements in the matrix needs to be chosen carefully
so to respect the same “wrap” boundary condition behaviour that we maintained in the primal
domain
1
.
1 The shifting property of convolutions can give you a hint of exactly how to structure and place the non-zero kernel
matrix elements in , before taking its DFT.
Given the appropriately constructed , and its DFT , we can obtain the DFT of the sampled output
image (i.e., with ) as
where is an element-wise product. Taking the real component of the inverse DFT of , we arrive at the
final output blurred image matrix .
𝐼(𝑥, 𝑦) 𝑘(𝑥, 𝑦)
𝐵(𝑥, 𝑦) = 𝐼(𝑥, 𝑦) ⊛ 𝑘(𝑥, 𝑦) 𝐼(𝑥, 𝑦) ≡ 𝐼[𝑥, 𝑦]
𝑘(𝑥, 𝑦) ≡ 𝑘[𝑥, 𝑦]
𝐵(𝑥, 𝑦) ≡ 𝐵[𝑥, 𝑦]
𝐈 (𝐈)𝑖,𝑗 = 𝐼(𝑖,𝑗) ˆ𝐈
5 𝐊 ˆ𝐊
𝐊 𝐈
𝐊
Deliverable 5 [25%]
Fourier Kernel Matrix: Complete the implementation of the FourierKernelMatrix that takes as
input an input blur kernel (with odd ), a 2-tuple of the shape (i.e., dimensions) of the
input/output image (image_shape), and (optionally) a forward 2D DFT function (DFT2D_func)
that you will call in your FourierKernelMatrix solution. The function returns the matrix
corresponding to the DFT of . The non-zero kernel elements in should be placed in the
matrix so to admit the same post-convolution wrap-based blurring behaviour as we saw in
Assignments 1 and 3.
☑
𝑛 × 𝑛 𝑛
ˆ𝐊
𝐊 𝐊
When implementing and debugging FourierKernelMatrix, you don’t necessarily have to rely on
your DFT2D routine (e.g., by passing something other than DFT2D as the DFT2D_func parameter).
ⓘ
𝐊
𝐊 ˆ𝐊 ˆ𝐁
𝐁 (𝐁)𝑖,𝑗 = 𝐵(𝑖,𝑗)
ˆ𝐁 = ˆ𝐈 ⋆ ˆ𝐊 ,
⋆ ˆ𝐁
𝐁
4.2 2D Fourier Deconvolution
In an unregularized setting, we can readily express the deconvolution problem using the DFT as
where is an element-wise division. Taking the real component of the inverse DFT of yields the
unblurred image matrix . Since the conditioning of the DFT and iDFT are 1, the underlying conditioning of
the convolution/deconvolution system remain as (un)stable as they were in the primal domain.
4.2.1 Laplacian Regularization
We provide you with a blurred image polluted by mild noise, as well as an blur kernel (which you will
need to extend appropriately into , as above).
As mentioned briefly in Assignment 3, Laplacian regularization is well-suited to image deblurring problems:
it adds a regularization term that penalizes large second derivatives in image-space, a statistical property of
many natural images.
In the primal domain, the Laplacian regularization kernel is
which you can extend and arrange into a matrix with dimensions that match the image resolution and with
structure that respects the wrapping boundary condition treatment (i.e., similar to the blur kernel’s matrix
structure).
We can solve the regularized least squares deconvolution problem, expressed in the primal domain,
in the frequency domain as
ˆ𝐈 = ˆ𝐁 ⋆ ,
−1 ˆ𝐊 (6)
⋆
−1 ˆ𝐈
𝐈
Unlike Assignments 1 and 3, we will not be providing you with the ground truth deblurred image for
your reference, below. You can code up example tests for forward convolution and nonnoisy/unregularized deblurring to sanity check your DFT/iDFT implementations; your code from A3 can
come in handy, here!
𝑛 × 𝑛
𝐊
3 × 3
𝑙(𝑥, 𝑦) ≡ 𝑙[𝑥, 𝑦] =
⎡
⎣
⎢
⎢
0
1
0
1
−4
1
0
1
0
⎤
⎦
⎥
⎥
(7)
𝐋
𝐊
min ‖𝐵(𝑥, 𝑦) − 𝐼(𝑥, 𝑦) ⊛ 𝑘(𝑥, 𝑦) + 𝜆‖𝐼(𝑥, 𝑦) ⊛ 𝑙(𝑥, 𝑦)
𝐼
‖
2
2 ‖
2
2
(8)
where is the DFT of and the absolute value of a complex value is . Taking the
real component of the inverse DFT of yields our regularized and deblurred image matrix . Note that in
Equation the regularized kernel is already in an inverted form, which is why an element-wise
multiplication ( ) is employed instead of an element-wise division ( , as in Equation ); i.e., with ,
Equations and are equivalent.
5 You’re Done!
Congratulations, you’ve completed the 4
th assignment. Review the submission procedures and guidelines
before submitting your Python script file assignment solution. Recall that any test code you incude in your
mainline (__main__) will not be graded, however it must still respect the collaboration/plagiarism polices.
formatted by Markdeep 1.13
ˆ𝐈 = ˆ𝐁 ⋆ ,
⎛
⎝
⎜
⎜
⎜
ˆ𝐊
+ 𝜆
∣
∣ˆ𝐊∣
∣
2
∣
∣ˆ𝐋∣
∣
2
⎞
⎠
⎟
⎟
⎟
ˆ𝐊
−1
reg
(9)
ˆ𝐋 𝐋 | ⋅ | 𝑎 + 𝑏ı 𝑎
2 + 𝑏 √‾‾‾‾‾‾‾2
ˆ𝐈 𝐈
9 ˆ𝐊
−1
reg
⋆ ⋆
−1 6 𝜆 = 0
6 9
Deliverable 6 [30%]
Laplacian-regularized Fourier Deconvolution Kernel: Complete the implementation of
LaplacianInverseFourierKernel to form and return . The function parameters match
those of FourierKernelMatrix. Choose and hardcode a regularization coefficient value that
yields a qualitatively good deblurred image (i.e., when applied to our blurred_noisy_image test
data using Equation ; see our test code in __main__: part of your grade for this deliverable will
be based on whether we can discern important details from the deblurred image.
☑
ˆ𝐊
−1
reg
𝜆
9
✒
