Update ST2D_exercise_withSolutions.ipynb

0310f4c9 · monj · c6bc7e15 · 0310f4c9
Commit 0310f4c9 authored 3 years ago by monj
--- a/ST2D_exercise_withSolutions.ipynb
+++ b/ST2D_exercise_withSolutions.ipynb
@@ -401,7 +401,7 @@
    }
   ],
   "source": [
-    "filename = './example_data_2D/OCT_im_org.png';\n",
+    "filename = './example_data_2D/oct_im_org.png';\n",
    "sigma = 0.5\n",
    "\n",
    "interact(st_and_hists, sigma = fixed(sigma), rho = (0,10,1), filename = fixed(filename));"

 %% 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 this exercise, you will be asked to write some code here and there, please write the code between the comments **###USER INPUT and ###END OF USER INPUT**. Note that the code is already written in the "_withSolutions" version of the exercise.
 Run the script cell by cell by placing the cursor inside the cell and pressing Run or Shift+Return.
 %% 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.
 - [structure_tensor](https://github.com/Skielex/structure-tensor) contains the structure tensor implementation that you will install using the *environment.yml* file.
 - *utilsST_2D* is a Python file containing helper functions for this exercise
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import skimage.io
 import matplotlib.pyplot as plt
 from structure_tensor import eig_special_2d, structure_tensor_2d
 from utilsST_2D import *
 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 <br>
 `S = structure_tensor_2d(img, sigma, rho)`
 - Compute the eigen values for the principal orientations of your structure, and the dominant orientation vector (**vec** $= \text{vec}[0]\hat{x}+\text{vec}[1]\hat{y}$). <br>
 `val,vec = eig_special_2d(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: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 `structure_tensor_2d` and `eig_special_2d` functions described above
 2. Calculate the orientation angles
 > - Compute the angles[rad] from the components of the dominant orientation (__vec__). <br>
 >> Look at the figure above and use a trigonometric function: `np.arccos()`, `np.arcsin()` or `np.arctan2()`.
 %% Cell type:code id: tags:
 ``` python
 ## RUN THE STRUCTURE TENSOR
 # Read the image
 filename = './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 = structure_tensor_2d(image.astype('float'), sigma, rho) # compute structure tensor matrix S
 val, vec = eig_special_2d(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]) # compute angles
 ### END USER INPUT ###
 # Prepare the angles to show as image
 angle /= np.pi #scale from 0 to 1
 orientation_image = angle.reshape(image.shape) # shape as an image
 # 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, 50, 50) array.
    Orientation information is carried in a (2, 50, 50) 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')
 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**. In this way, we will keep the orientation angle for highly anisotropic structures (anisotropy of a line = 1), downweigh the angle as it becomes less reliable, and not show it at all for a circular structure (anisotropy of a circle = 0).
 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(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 = './example_data_2D/10X.png';
 sigma = 0.5
 interact(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 = './example_data_2D/OCT_im_org.png';
+filename = './example_data_2D/oct_im_org.png';
 sigma = 0.5
 interact(st_and_hists, sigma = fixed(sigma), rho = (0,10,1), filename = fixed(filename));
 ```
 %% Output
 %% 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)