# Lecture 13:  Integral Transforms, D/FFT and Electron Microscopy
----

## Reading and Reference

* Advanced engineering Mathematics, E. Kreyszig, John wiley and Sons, 2010
* Numerical Recipes, W. Press, Cambridge University Press, 1986
* M. De Graef and M. McHenry, Structure of Materials, Cambridge University Press, 2nd ed.
* C. Hammond, The Basics of Crystallography and Diffraction, Oxford Science Publications, 4th ed.

## What to Learn?

* The definition of an integral transform
* The algorithm for computing the discrete Fourier transform
* How diffraction patterns can be used to create phase contrast images in electron microscopy

## What to Do?

* Compute the Fourier transform of different aperture functions.
* Select different regions of a Fourier transform to reconstruct a simulated TEM image and an image of your choosing.

### Introduction to Integral Transforms
----

An integral transform maps a function of one independent variable into a function of another independant variable using a _kernel_:

$$g(\alpha) = \int_{a}^{b} f(t) K(\alpha,t) dt $$

The function $f(t)$ is transformed to a new function $g(\alpha)$ through the definate integral.  A similarity to the dot product of functions is evident in this form and this operation can be thought of as a mapping or projection of $f(t)$ into a different independent variable $\alpha$.  Existence, integrability and inversion of integral transform operations are important in the study of this topic, although not covered in these notes.

Two examples of integral transforms, the Laplace and Fourier, are discussed in this lecture.  It is typical to use the Laplace transform to remove the time dependence from Fick's second law in diffusion problems.  The Fourier transform is used in the study of diffraction under certain conditions.

To assist in this lecture some special symbols in `Python` and `sympy` are reviewed:


```python
import sympy as sp
sp.init_printing(use_latex=True)
```


```python
# symbols we will need below
x,y,z,t,c = sp.symbols('x y z t c')
# note the special declaration that omega is a positive number
omega = sp.symbols('omega', positive=True)
```

###  Complex Number Review
----

A reminder that $i$ is the square root of negative one and this is how you specify $i$ in `Sympy` and that is different than the complex data type in `Python`.


```python
sp.I**2
```

The natural logarithm of $e$ is $1$:


```python
sp.log(sp.E)
```

In SymPy there are two ways to deal with integration.  If you would like to represent an unevaluated integral, you can use the `Integral` function.  If you want to compute the integration of an expression you can use the `integrate` function.


```python
sp.Integral(sp.E**(sp.I*omega*t),t)
```


```python
# 'omega', positive=True
sp.integrate(sp.E**(sp.I*omega*t),t)
```

Where we assume there is no zero frequency (as we are dividing by $\omega$) - hence the assumption `positive=True` in the symbol definition above.  (Try replacing $\omega$ with $y$ and inspect the expression returned by `integrate`.)

### The Fourier Transform
----

As the domain of the periodicity increases, the frequency spectrum required required to represent the function becomes more finely divided.  Recall the argument of the trigonometric terms in the functions of the Fourier series:

$$ \frac{n \pi (\omega +c)}{d} $$

where n is the order of the frequency component, c the offset relative to the origin, and d the domain width.  If we let the domain width go to infinity (implying that the function is not periodic) then an integral sum is required rather than a discrete summation.  The, infinte, non-periodic function and its frequency spectrum are related by the Fourier transform defined by:

$$ \hat{f}(\omega) = \sqrt{\frac{1}{2\pi}} \int^{+\infty}_{-\infty} f(t) \exp[-i \omega t] dt $$

This results in a mapping of the function f(t) into frequency space.  

The real or complex and even or odd nature of the function $f(t)$ determines if the transformed function is even, odd, real, or complex.  For the purposes of materials crystal structures in this lecture we will be using even and real functions.  

### Diffraction from An Aperture
----

A useful physical problem requiring use of the Fourier transform is diffraction.  In this problem we will use a top-hat function to represent the location of an infinity of wave sources from an aperture.  We use the `sp.Piecewise` function to generate a "tophat" function for the Fourier transform.


```python
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 8, 4

p = sp.Piecewise((0,x<-1),(1,x<1),(0,True))
sp.plot(p);
```

At some distance from the aperture we place a detector that measures the combined intensity of all the wave sources, however due to the finite width of the slit each wave travels a different distance to the detector.  The phase difference between the waves at the detector is given by the Fourier transform of the aperture function when the [Fraunhofer](https://en.wikipedia.org/wiki/Fraunhofer_diffraction_equation) approximation is valid.

This aperture function is even and real so we expect our transformed function to also be even and real.  We use the definition of the integral transform above to write an explicit integral statement of the Fourier transform of the top-hat function above.  The integral is $1$ between $c$ and $-c$ and zero elsewhere - so we can integrate **just the non-zero part**.  This is integrated as: 


```python
sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,-c,c))
```

Calling explicitly for the integration and assigning the result to `a`:


```python
a = sp.sqrt(1/(2*sp.pi))*sp.integrate(1*sp.exp(-sp.I*2*omega*x),(x,-c,c))
a
```

This does not (at first glance) appear to be a real function due to the two exponential terms, but we can use some of the algebraic methods built into `SymPy` to help.  We can ask for this form using sines and cosines with the `rewrite` method.  Furthermore - we can simplify it further with the expand function.  Trial and error may be required to determine the best combination and ordering of algebraic manipulations.


```python
solution = sp.expand(a.rewrite(sp.sin))
solution
```

Here we can use the `subs` (substitution) method to set the value of `c`.  I plotted the square of the function since the intensity of a diffracted wave is related to the time averaged energy transferred by the wave.  This is proportional to the amplitude squared.  As our function is real valued, we can just plot the square.


```python
sp.plot(solution.subs(c,1));
```


```python
sp.plot(solution.subs(c,1)**2);
```

### Diffraction from Two Apertures
----

We could perform the same integration over two top-hat functions and plot those results.


```python
compositeIntegral = sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,1,2)) + \
sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,-2,-1))
compositeIntegral
```


```python
om = compositeIntegral.doit()
om
```

The diffracted intensity from this pair of slits would appear as:


```python
sp.plot(om.rewrite(sp.sin).expand()**2)
```

Or we could functionalize our function to explore other parameters:


```python
def diffractionFunction(d=4.0, w=1.0):
    result = sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),\
                                     (x,-(d+w),-(d-w))) + \
    sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),\
                                     (x,(d-w),(d+w)))
    return result.doit()
```


```python
sp.expand(diffractionFunction(10.,2.).rewrite(sp.sin))
```

### DIY:  Complex Numbers
----

Perform the Fourier transformation on an odd or complex valued function.  Plot the real and imaginary parts of both the target function and the transformed functions.

### DIY:  The Airy Disk
----

Solve for the diffracted intensity in two dimensions from a circular aperture. It may be easier to do this as a discrete problem using the DFT below. 

### The Discrete Fourier Transform
----

The discrete Fourier Transform is defined [here](http://en.wikipedia.org/wiki/Discrete_Fourier_transform) and is regarded as one of the most important advances in computing science in the 20th century.  Other resources such as Numerical Recipes, the Python help files and many other websites detail the calculation and implementations.

It is often instructive to review other implementations of the DFT to help you gain experience.  I will be modeling this implementation after Jake Vanderplas' blog article [here](http://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/).  Following the notion in the blog article:

Forward DFT:

$$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}$$

Inverse DFT:

$$x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N}$$

In this section of the notebook, we use Vanderplas' description and implementation.

----
For simplicity, we'll concern ourself only with the forward transform, as the inverse transform can be implemented in a very similar manner.  Taking a look at the DFT expression above, we see that it is nothing more than a straightforward linear operation: a matrix-vector multiplication of $\vec{x}$,

$$\vec{X} = M \cdot \vec{x}$$

with the matrix $M$ given by

$$M_{kn} = e^{-i~2\pi~k~n~/~N}$$

With this in mind, we can compute the DFT using simple matrix multiplication as follows:


```python
import numpy as np
def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)
```

We can use the "all close" function to check if the result from `DFT_slow` and `Numpy` are close:


```python
x_signal = np.random.random(1024)
np.allclose(DFT_slow(x_signal), np.fft.fft(x_signal))
```

----
 

I think it would be instructive to symbolically expand the matrix above so that it is clear how `n*k` leads to a two dimensional matrix.  Switching to `sympy` symbols to expose the details we can do the following:


```python
import sympy as sp
from sympy import Matrix
import numpy as np
sp.init_printing()
```

* `x` is the input vector.
* `k` is the wavenumber or frequency.
* `n` is the component of the input vector.


```python
x = sp.Matrix(sp.symbols('x0:5'))
n = sp.Matrix(sp.symbols('n0:5')).T
k = sp.Matrix(sp.symbols('k0:5'))
N = sp.symbols('N')
M = (-sp.I*2*sp.pi*k*n/N).applyfunc(sp.exp)
```


```python
M*x
```

Each frequency element is projected into each point of the input vector - the matrix links `k` and `n`.  So - the contribution at each point is a sum of each frequency contribution, similar to the dot product of functions.

### DFT with Numpy Functions
----

In this section we use the `FFT` submodule of `numpy` to help in the computation of the DFT.


```python
?np.fft # This gives us information on the conventions used in the return values of the functions.
```


```python
?np.fft.fft # This is the main DFT function we will use.
```


```python
?np.fft.fftfreq  # This is a helper function to prepare a vector of frequencies.
```


```python
?np.arange  # Points in an evenly spaced interval.
```

This approach is derived from a nice discussion on FFT found on the blog Glowing Python.

First we will divide up time into `samplingInterval` sized chunks between 0 and 1.  This will aid in getting the x-axis scaled correctly so that frequency can be read directly off the DFT result.  You can take `samplingInterval` in seconds putting samplingRate in Hz.  Notice the approach here - we could have done this all in one line, but, by intelligently naming our variables and exposing the details of our thoughts the code is more readable:


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

samplingRate = 150.0
samplingInterval = 1.0/samplingRate
timeVector = np.arange(0, 1, samplingInterval)

# Print out the first few elements so you can see what is going on:
timeVector[0:10:]
```

Next we decide on the frequency of our signal and create a list to have a signal to work with. 


```python
signalFrequency = 10.0;
ourSignal = np.sin(2*np.pi*signalFrequency*timeVector) + 0.5*np.sin(2*np.pi*(2*signalFrequency)*timeVector)
```

Plotting the input function for clarity:


```python
fig = plt.figure()

axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
axes.plot(timeVector, ourSignal, 'r')
axes.set_xlabel('Time')
axes.set_ylabel('Signal')
axes.set_title('Our Modest Signal');
```

Using `numpy` to compute the DFT:


```python
n = ourSignal.size
frequencies = np.fft.fftfreq(n, d=1.0/samplingRate)
spectrum = np.abs(np.fft.fft(ourSignal))

fig = plt.figure()

axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])

axes.scatter(frequencies, spectrum, c='r', marker='s', alpha=0.4)

axes.set_xlabel('Frequency')
axes.set_ylabel('Amplitude')
axes.set_title('Our Amplitude Spectrum');
```

### Interactive Microscopy Demonstration (Optional)
Original developed by C. Carter, translated to Python by D. Lewis

---

Transmission electron microscopy utilizes diffraction to determine crystal structures and develop contrast in images.  In this section of the lecture we will simulate the diffraction pattern of an atomic structure.  Using this diffraction pattern we will simulate using a diffraction aperture to reconstruct a phase contrast image.


```python
%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy.fft import *

def atomic_func(x,y):
    param = 64.0
    return (1+np.sin(4*(x+y)*2*np.pi/param))*(1+np.sin(2*(x-2*y)*2*np.pi/param))/4

def aperture(X, Y, xoffset, yoffset, size):
    return (X-xoffset)**2+(Y-yoffset)**2 > size**2
```

We define two functions above:

* `atomic_func` is used to provide an image function periodic in two dimensions from which the diffraction pattern will be constructed.  This can be thought of as the density of electrons in a solid that is used to approximate a crystal structure.
* `aperture` returns a Boolean array that will be used to mask the diffraction pattern so that individual frequencies can be selected for image reconstruction.  `aperture` will return `True` or `False`.


```python
x = np.arange(0.0,256.0,1.0)
y = np.arange(0.0,256.0,1.0)
X,Y = np.meshgrid(x, y)

Z = atomic_func(X,Y)
```

The `Z` array holds the atomic image function.


```python
P = np.zeros(Z.shape,dtype=complex)
K = np.zeros(Z.shape,dtype=complex)

K = fftshift(fft2(Z, norm='ortho'))

P = np.copy(K)
P[np.where(aperture(X, Y, 128, 128, 3) & aperture(X, Y, 150, 128, 3))] = 0
```

The `P` array holds the processed Fourier spectrum.  The values of `P` are set to zero when they are outside the aperture.  We use the `K` array to hold a opy of the image


In this cell we create two more `numpy` arrays (there are other ways to do this) that have the same shape as Z. The `P` array we use to hold the processed Fourier spectrum.  The processing uses `numpy`'s Boolean indexing to set values in P equal to zero when they are "outside" the aperture.  When we get to the images below you'll see what is meant.

Because Python passes by reference we need to call for a copy of K so that we can modify one without changing the other.

From this processed spectrum we will create an image.  The K array holds the whole Fourier spectrum.


```python
Im = fftshift(ifft2(P))
```

Above we reprocess `P` into the image `Im`.


```python
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(30,9))
axes[0].imshow(Z, origin='lower')
axes[1].imshow(abs(K),origin='lower', cmap=plt.get_cmap('pink'))
aperture1 = plt.Circle((128,128),3**2,color='r', fill = False)
aperture2 = plt.Circle((150,128),3**2,color='y', fill = False)
axes[1].add_artist(aperture1)
axes[1].add_artist(aperture2)
axes[2].imshow(abs(Im)**2, origin='lower')
plt.show()
```

### Homework
----

Apply the DFT to an image of your choosing.  Select the low frequency part of the DFT and regenerate the image (i.e. take the inverse FFT) from only these selected frequencies.  Use a Boolean selection to zero out parts of the frequency spectrum before you convert back.  To read an image in from disk, use the `ndimage` function from SciPy:

```python
from scipy.ndimage import imread
img = imread('./images/pattern2.jpg', mode='L')
```

checking the data type of `img` will prove helpful.

### Summary
----

* Integral transforms map one function space into another function space.  You can find books that include tables of Laplace and Fourier transforms.  Many other transforms exist - but the principle is the same.
* The DFT organizes amplitude information in predictable yet non-intuitive ways.  Read the documentation for the functions you use!  
* Integral transforms are a means for reducing the complexity of certain ODEs and PDEs.
* Diffraction and diffusion are two example applications where integral transforms can be employed.

### Reading Assignments and Practice
----

* Pam Champness' book on electron diffraction is a (relatively) easy read on diffraction.  You can always have a look at Cullity, Hammond, or any other book on structure and X-ray/electron characterization.
* Practice taking the FFT of signals you construct by hand.  This is a good step when you are debugging a problem.  You should always have a test case available to determine if your work is correct.

b133e5e87cfe8e5c42ee68ae9cac213947fba2e9
ipynb
Jupyter Notebook
Lecture-13-Integral-Transforms.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Lecture-13-Integral-Transforms.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Lecture-13-Integral-Transforms.ipynb
juhimgupta/MTLE-4720
41797715111636067dd4e2b305a782835c05619f
MIT
Qwen/Qwen-72B
1. YES
2. YES