Skip to content
Snippets Groups Projects
Commit 7269d1af authored by monj's avatar monj
Browse files

Updated structure tensor exercises

parent 9c59b184
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Structure Tensor/StructureTensor2D_exercise.ipynb
+ 200
− 81
View file @ 7269d1af
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Structure Tensor/StructureTensor2D_exercise_solutions.ipynb
+ 51
− 33
View file @ 7269d1af
%% Cell type:markdown id: tags:
# 2D Structure tensor - workshop exercise
The aim of this exercise is to get you familiarised with the structure tensor. You are going to learn how to:
- Use the functionality of the structure tensor.
- Obtain the orientation angle and degree of anisotropy from the structure tensor.
- Visualise the output of the structure tensor.
- Tune the input parameters of the structure tensor.
In some parts of code you will be asked to give user input, write in between the comments ###USER INPUT and ###END OF USER INPUT.
%% Cell type:markdown id: tags:
## Packages
Let's import the packages we will use in this exercise. The most important ones for you to know are:
- [numpy](www.numpy.org) is the basic package for scientific computing with Python.
- [skimage](https://scikit-image.org/) is one of the top image processing libraries in Python. Today, we will use its intput/output module io to read in the images, and the transform module to...
- [matplotlib](http://matplotlib.org) is the most commonly used package for plotting in Python.
- [st2d](https://lab.compute.dtu.dk/QIM/structure-tensor) contains the structure tensor implementation for 2D images by V.A. Dahl (vand@dtu.dk).
- *utilsST* is a Python file containing helper functions for this exercise
%% Cell type:markdown id: tags:
First, run the cell right below to get the structure tensor repository.
If this doesn't work, ignore the cell and instead do the following.
1. Download the [structure tensor repository](https://lab.compute.dtu.dk/QIM/structure-tensor) as a .zip.
2. Extract and rename to 'structure_tensor'
3. Place folder in the directory of this exercise.
%% Cell type:code id: tags:
``` python
!git clone "https://lab.compute.dtu.dk/QIM/structure-tensor"
import os
os.rename('structure-tensor','structure_tensor')
```
%% Output
Cloning into 'structure-tensor'...
remote: Enumerating objects: 96, done.
remote: Counting objects: 100% (96/96), done.
remote: Compressing objects: 100% (59/59), done.
remote: Total 96 (delta 33), reused 89 (delta 28), pack-reused 0
Unpacking objects: 100% (96/96), done.
Checking connectivity... done.
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
<ipython-input-3-28c94fc19bd1> in <module>
1 get_ipython().system('git clone "https://lab.compute.dtu.dk/QIM/structure-tensor"')
2 import os
----> 3 os.rename('structure-tensor','structure_tensor')
4 from structure_tensor import st2d
OSError: [Errno 66] Directory not empty: 'structure-tensor' -> 'structure_tensor'
%% Cell type:code id: tags:
``` python
import numpy as np
import skimage.io
#import scipy.ndimage
#import skimage.transform
import matplotlib.pyplot as plt
from structure_tensor import st2d
import utilsST
from ipywidgets import interact, IntSlider, fixed
```
%% Cell type:markdown id: tags:
## Structure tensor functions
The main functionality of the structure tensor is implemented in the two functions below:
- Compute the Structure Tensor for 2D Data
S = st2d.structure_tensor(img, sigma, rho)
- Compute the eigen values for the principal orientations of your structure, and the dominant orientation vector ($vec[0]\hat{x}+vec[1]\hat{y}$).
val,vec = st2d.eig_special(S)
<img src="structuretensor.png" width="100%">
*To get HELP ON A FUNCTION, write funcion_name? in an empty cell.*
*To INSERT AN EMPTY CELL, press the plus symbol +, second icon on the left underneath File.*
%% Cell type:code id: tags:
``` python
st2d.structure_tensor?
```
%% Cell type:markdown id: tags:
## Part 1. Guided implementation of the structure tensor
This first part will guide you through the implementation of your first analysis using structure tensor. You are going to write code to:
1. Compute the structure tensor matrix and the orientation vectors.
Use the st2d functions described above
2. Calculate the orientation angles and reshape back to an image.
- Compute the angles[°] from the dominant orientations (vec).
Use np.pi and one of the trigonometric functions: np.arccos(), np.arcsin() or np.arctan2(). Look at the figure above.
- Reshape your output with output.reshape(img.shape)
%% Cell type:code id: tags:
``` python
## RUN THE STRUCTURE TENSOR
# Read the image
filename = './structure_tensor/example_data_2D/drawn_fibres_B.png'; # path of your image
image = skimage.io.imread(filename) # read
print(f'The image has a shape {image.shape}, i.e. {image.size} pixels.') # print image dimensions
# Initialise parameters
sigma = 0.5
rho = 2
# Compute the structure tensor
### USER INPUT HERE ### (≈ 2 lines of code)
S = st2d.structure_tensor(image, sigma, rho) # compute structure tensor matrix S
val, vec = st2d.eig_special(S) # compute dominant orientation vector (vec) and principal orientation weights (val)
### END USER INPUT ###
print(f'Structure tensor information is carried in a {S.shape} array.')
print(f'Orientation information is carried in a {vec.shape} array.')
# Compute orientation angles
### USER INPUT HERE ### (≈ 1 line of code)
angle = np.arctan2(vec[1], vec[0])/np.pi # compute angles
orientation_image = angle.reshape(image.shape) # shape as an image
### END USER INPUT ###
# Visualise the results
figsize = (10,5)
fig, ax = plt.subplots(1, 2, figsize=figsize, sharex=True, sharey=True)
ax[0].set_title('Input image')
ax[0].imshow(image, cmap='gray')
ax[1].set_title('Orientation as color')
ax[1].imshow(orientation_image, cmap = 'hsv')
plt.show()
```
%% Output
The image has a shape (50, 50), i.e. 2500 pixels.
Structure tensor information is carried in a (3, 2500) array.
Orientation information is carried in a (2, 2500) array.
%% Cell type:markdown id: tags:
With the structure tensor, we get an orientation value for every pixel, but actually, we are only interested in the orientations of the fibres.
Because the fibres in this image are bright and the background is dark, the orientations can be overlaid on the image by multiplying the image pixels by the orientation values (see below).
We will also use the function `st2d.plot_orientations()` to the visualise the orientations as vectors.
%% Cell type:code id: tags:
``` python
## VISUALISE THE ORIENTATIONS OVERLAID ON THE IMAGE
figsize = (10,5)
fig, ax = plt.subplots(1, 2, figsize=figsize, sharex=True, sharey=True)
ax[0].set_title('Orientation as arrows on image')
ax[0].imshow(image, cmap='gray')
st2d.plot_orientations(ax[0], image.shape, vec)
ax[1].set_title('Orientation as color on image')
ax[1].imshow(plt.cm.gray(image)*plt.cm.hsv(orientation_image))
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
When visualizing only the dominant orientation, it may seem as the orientation changes abruptly in some areas (e.g. where the fibres merge). These areas where the orientation is smoothly changing from one dominant orientation to another dominant orientation, are actually isotropic, meaning that they do not actually have a preferencial orientation.
%% Cell type:markdown id: tags:
## Part 2. Shape of the Structures
While some image regions may have a preferencial orientation, others don't. In this part, we are going to measure the isotropy of the structures to find out whether they have a preferencial direction. Are they more linear-like (isotropy of a line = 0) or circular-like (isotropy of a circle = 1)? We will use this information to give more importance to the angles in regions with preferential orientations.
The first output of the structure tensor (`val`) indicates how fast the intensities change around a pixel in the two principal directions. Along the first principal direction of a line, the intensities will not change, and `val[0]` will be very small. Along the second principal direction, the intensities will change fast so `val[1]` will be large. For a circle, both values will be equal.
The ratio between the components of `val` will give us the **degree of structural isotropy**. To obtain a more realistic idea of the orientations of our material, we will weigh the dominant orientation by the anisotropy.
Can you find the formula for the anisotropy, given `val[0]` and `val[1]`?
%% Cell type:code id: tags:
``` python
## WEIGH THE ORIENTATIONS BY THE ANISOTROPY
### USER INPUT
anisotropy = 1-val[0]/val[1] # compute anisotropy
### END OF USER INPUT
anisotropy = anisotropy.reshape(image.shape)
# Plot image, anisotropy and weighed orientations
fig,ax = plt.subplots(2, 3, figsize=(15,10), sharex=True, sharey=True)
ax[0][0].imshow(plt.cm.gray(image))
ax[0][0].set_title('Image')
ax[0][1].imshow(anisotropy)
ax[0][1].set_title('Degree of anisotropy')
ax[0][2].imshow(plt.cm.gray(anisotropy)*plt.cm.hsv(orientation_image))
ax[0][2].set_title('Orientation and anisotropy')
# Overlaid on the images
ax[1][0].imshow(plt.cm.gray(image)*plt.cm.hsv(orientation_image))
ax[1][0].set_title('Orientation on image')
ax[1][1].imshow(image*anisotropy)
ax[1][1].set_title('Degree of anisotropy on image')
ax[1][2].imshow(plt.cm.gray(image)*plt.cm.gray(anisotropy)*plt.cm.hsv(orientation_image))
ax[1][2].set_title('Orientation and anisotropy on image')
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
## Part 3. Understand the parameters sigma and rho
The aim of this part of the exercise is to get you familiarised with the parameters sigma and rho.
In the structure tensor, there is a first step to remove noise and a second step to compute the orientations, where the dominant direction is found as the direction in which the intensities change the least. In the first step, noise removal, the image is smoothed with a kernel of size sigma. In the second step, the intensities are integrated around the pixel in a region of size defined by rho.
Now that you know that sigma is related to the scale of the image noise and rho to the scale of the structures of interest, play with the widgets to change the values of the parameters and see the effects on the results. Move the widgets and wait for the image to update.
%% Cell type:code id: tags:
``` python
interact(utilsST.st_and_plot, sigma = (0.25,5,0.25), rho = (0.1,15,0.5), image = fixed(image));
```
%% Output
%% Cell type:markdown id: tags:
Set rho low (< 1) and play with the noise level (sigma in 0-1). See how the noise is reduced.
Increase rho and watch the orientations become more smooth as the size of the integration window covers the width of the structures. As rho is increased (try rho > 6), the smaller structures start to disapear in the anisotropy image and the orientations colours start to blend across structures.
%% Cell type:markdown id: tags:
## Part 4. Application: Find whether your sample has a preferencial orientation
The structure tensor is very useful for finding out whether your material has a preferencial orientation. This has numerous applications, from manufacturing of composite materials to understanding blood flow in your body.
We are going to see two examples, one about cardboard fibres and the other about retina vasculature. You will see how the orientation information can be accumulated into a histogram of angles and used to determine the preferencial orientations of your sample, if any!
%% Cell type:markdown id: tags:
### Orientation analysis of carboard fibres.
Here, the orientation information is collected from the whole image. However, as shown above, the orientation information is reliable only in areas of high anisotropy. For this reason, and depending on problem at hand, it may be desirable to weight the orientation using the anisotropy (remove parts with low anisotropy) and/or the intensity (remove the background and keep only information of fibres).
%% Cell type:code id: tags:
``` python
filename = './structure_tensor_master/example_data_2D/10X.png';
filename = './structure_tensor/example_data_2D/10X.png';
sigma = 0.5
interact(utilsST.st_and_hists, sigma = fixed(sigma), rho = (2,20,3), filename = fixed(filename));
```
%% Output
%% Cell type:markdown id: tags:
Note: To choose the scale at which you wish to capture the orientations, you may want to first analyse a smaller region of interest so as to visualise the data close-up and determine a sigma, and especially a rho parameter, that provides the orientations you are interested in. See the difference between choosing a rho of 17 and a rho of 2.
%% Cell type:markdown id: tags:
### Orientation information from the OCT image of a retina
%% Cell type:markdown id: tags:
Play with the widget and find the optimal rho value for capturing the orientation of the thicker vessels.
%% Cell type:code id: tags:
``` python
filename = './structure_tensor_master/example_data_2D/OCT_im_org.png';
filename = './structure_tensor/example_data_2D/OCT_im_org.png';
sigma = 0.5
interact(utilsST.st_and_hists, sigma = fixed(sigma), rho = (0,10,1), filename = fixed(filename));
```
%% Cell type:markdown id: tags:
%% Output
Play with the widget and find the optimal rho value for capturing the orientation of the thicker vessels.
%% Cell type:markdown id: tags:
When working with large amounts of data, the computational time becomes crucial. In [Github](https://github.com/Skielex/structure-tensor/tree/master/structure_tensor) you can find an optimised version of the structure tensor developed at the Technical University of Denmark (DTU )by Niels Jeppesen (niejep@dtu.dk). This implementation provides the option of using the package cupy, instead of numpy, to run the computations on a Graphical Processing Unit (GPU).
When working with large amounts of data, the computational time becomes crucial. In [Github](https://github.com/Skielex/structure-tensor/tree/master/structure_tensor) you can find an optimised version of the structure tensor developed at the Technical University of Denmark (DTU) by Niels Jeppesen (niejep@dtu.dk). This implementation provides the option of using the package cupy, instead of numpy, to run the computations on a Graphical Processing Unit (GPU).
*Please note that the visualisation functions used today would have to be adapted to work with Niels' code, as the output of his functions is in a different format.*
%% Cell type:markdown id: tags:
This workshop exercise was made by Monica J. Emerson (monj@dtu.dk) based on the example scripts by Vedrana A. Dahl (vand@dtu.dk)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment