03 - Orientation Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lineament Analysis\n",
    "\n",
    "Let's take what we've learned about edge filters and start to do something more interesting...\n",
    "\n",
    "Many of us have spent time trying to interpret, digitize, and analyze orientations of linear features that are visible in data of some sort.  There are very good reasons why it's often done by hand, but sometimes it's possible to automate.\n",
    "\n",
    "Let's take a look at some aerial photography data from near Arches National Park in Utah, USA.  There are prominent linear features (joints/fractures, in this case) visible in the imagery.  Ideally, we want a rose diagram of them without manually digitizing each one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio as rio\n",
    "import scipy.ndimage\n",
    "\n",
    "from context import data\n",
    "from context import utils\n",
    "\n",
    "with rio.open(data.naip.lineaments, 'r') as src:\n",
    "    aerial_photo = src.read()\n",
    "\n",
    "# Rasterio has a \"bands-on-first-axis\" convention, matplotlib/etc has\n",
    "# a \"bands-on-last-axis\" convention. Use moveaxis to switch between.\n",
    "aerial_photo = np.moveaxis(aerial_photo, 0, -1)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(aerial_photo)\n",
    "ax.set(xticks=[], yticks=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this type of analysis, we often don't need color information.  It's easiest to leave it out.  We'll analyze grayscale data instead.  Let's use a simple average of RGB values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gray_aerial = aerial_photo.astype(float).mean(axis=-1)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(gray_aerial, cmap='gray')\n",
    "ax.set(xticks=[], yticks=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our lineaments are basically edges. Remember that we said gradient magnitude is a type of edge detector?  Let's go ahead and look at the raw gradient magnitude. We don't care about absolute values at all in this case, so let's use a Sobel filter.\n",
    "\n",
    "Last time we used `scipy.ndimage.generic_gradient_magnitude` and `scipy.ndimage.sobel`.  However, that's a bit verbose and we don't care about the separate X and Y components, so let's use scikit-image's Sobel filter method instead, which is a bit simpler.  We'll use the \"toggler\" again so that we can easily compare it to the original imagery."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage.filters\n",
    "\n",
    "im_grad_mag = skimage.filters.sobel(gray_aerial)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(aerial_photo)\n",
    "im = ax.imshow(im_grad_mag, cmap='gray_r', label='Sobel', vmin=0, vmax=50)\n",
    "ax.set(xticks=[], yticks=[])\n",
    "utils.Toggler(im).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oy... That's noisy... \n",
    "\n",
    "Thankfully, we just talked about gaussian gradient magnitude as a way of producing less noisy gradients.  Let's apply it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 3\n",
    "gauss_grad_mag = scipy.ndimage.gaussian_gradient_magnitude(gray_aerial, sigma)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(aerial_photo)\n",
    "im1 = ax.imshow(im_grad_mag, cmap='gray_r', label='Sobel', vmin=0, vmax=50)\n",
    "im2 = ax.imshow(gauss_grad_mag, cmap='gray_r', label='Gauss', vmin=0, vmax=15)\n",
    "ax.set(xticks=[], yticks=[])\n",
    "utils.Toggler(im1, im2).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, we have something now that we could think about using directly.  Let's try thresholding the gaussian gradient magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_thresh = gauss_grad_mag > 5\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(gray_aerial, cmap='gray')\n",
    "im = ax.imshow(np.ma.masked_equal(grad_thresh, 0), vmin=0,\n",
    "               interpolation='nearest')\n",
    "ax.set(xticks=[], yticks=[])\n",
    "\n",
    "def update(thresh):\n",
    "    grad_thresh = gauss_grad_mag > thresh\n",
    "    im.set_data(np.ma.masked_equal(grad_thresh, 0))\n",
    "\n",
    "utils.Slider(ax, 1, 10, update, start=3).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can fairly easily imagine skeletonizing the classification we just created to get the exact edges.  It's the same idea as what we used to get a nice line at the toe of slope around seamounts earlier.\n",
    "\n",
    "You might even imagine trying to do a bit better job of skeletonizing so that nearby \"ridges\" linked up instead of being separate features.\n",
    "\n",
    "That's the basic idea behind the [Canny filter](https://en.wikipedia.org/wiki/Canny_edge_detector). \n",
    "\n",
    "It's essentially skeletonizing a thresholding gaussian gradient magnitude, but it tries to join up nearby features to give nice, continuous lines.  Let's give it a try on this data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage.feature\n",
    "\n",
    "canny = skimage.feature.canny(gray_aerial, sigma=3)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(gray_aerial, cmap='gray')\n",
    "ax.imshow(np.ma.masked_equal(canny, 0), vmin=0, interpolation='nearest')\n",
    "ax.set(xticks=[], yticks=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay! Now we're getting closer to identifying our lineaments.  However, we've still got a lot of noise and a lot of small features we're not too interested in.\n",
    "\n",
    "Not only that, but how are we going to turn these into actual lineaments we can get an orientation of?  If we vectorize this result, the line segments would be mostly curved and not pointing along the features we're interested in.\n",
    "\n",
    "Basically we want to try to extract straight lines from these ridges.\n",
    "\n",
    "A good way to detect straight lines in an image is to use a [Hough Transform](https://en.wikipedia.org/wiki/Hough_transform).  The basic idea is to progressively rotate the image and sum along rows of the rotated result, then add each summed version as a new column.  Straight lines will form a local peak in the resulting array.  I'm not going to go over this in too much detail, though.\n",
    "\n",
    "The key thing to know about a Hough transform is that it gives an _infinite_ line. I.e. `y = Ax + B`.  You don't get start and end points, which is usually what we're interested in.\n",
    "\n",
    "Therefore, there's a variant called a [probabilistic Hough transform](https://scikit-image.org/docs/dev/auto_examples/edges/plot_line_hough_transform.html#probabilistic-hough-transform) that attempts to identify likely straight line segments of a given length.  It results in discrete line segments with specific locations.  The length can be used as \"knob\" to tune whether you're finding lots of small linear features or fewer large linear features.  \n",
    "\n",
    "Sounds perfect for this task! Let's appy it to our Canny-filtered edges above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import probabilistic_hough_line\n",
    "from matplotlib.collections import LineCollection\n",
    "\n",
    "gap_ratio = 0.12\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(aerial_photo)\n",
    "col = ax.add_collection(LineCollection([], color='yellow'))\n",
    "ax.set(xticks=[], yticks=[])\n",
    "\n",
    "def update(length):\n",
    "    lines = probabilistic_hough_line(canny, line_length=length, \n",
    "                                     line_gap=int(length*gap_ratio))\n",
    "    col.set_segments(lines)\n",
    "\n",
    "utils.Slider(ax, 5, 50, update, start=10).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No matter how we do this, we pick up some features we're not interested in, and leave out some that we are.  Overall, it's reasonable if we're mostly interested in orientations.  Let's show a rose diagram of the lineaments we identified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "lines = probabilistic_hough_line(canny, line_length=30, \n",
    "                                 line_gap=5)\n",
    " \n",
    "# Calculate azimuth\n",
    "lines = np.array(lines)\n",
    "dx, dy = np.squeeze(np.diff(lines, axis=1)).T\n",
    "# Negative dy due to image orientation, 90 - angle for azimuth\n",
    "angles = np.pi / 2 - np.arctan2(-dy, dx)\n",
    "\n",
    "fig = plt.figure(constrained_layout=True)\n",
    "ax1 = fig.add_subplot(2, 1, 1)\n",
    "ax2 = fig.add_subplot(2, 1, 2, projection='polar', theta_offset=np.pi/2,\n",
    "                      theta_direction=-1)\n",
    "ax1.imshow(aerial_photo)\n",
    "ax1.add_collection(LineCollection(lines, color='yellow'))\n",
    "ax2.hist(np.concatenate([angles, angles + np.pi]), bins=60)\n",
    "\n",
    "ax1.set(xticks=[], yticks=[])\n",
    "ax2.set(xticks=[], yticks=[], axisbelow=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, this worked. However it did a poor job of identifying many of the prominent lineaments.  It also identifies a lot of lineaments that are purely perpendicular to the sun direction. We should be able to do better.\n",
    "\n",
    "One of the ways we can identify lineaments is to look for _consistent directions_ of image gradients.  Rather than just looking at the magnitudes, let's take direction into account as well.\n",
    "\n",
    "A common and very useful technique is the [structure tensor](https://en.wikipedia.org/wiki/Structure_tensor). The idea should be familiar to most folks who've worked on orientation statistics.  We take the gradient vectors of the image within a moving window for each pixel.  Then we build a 2x2 covariance matrix from the dx, dy componenents of those gradient vectors.  This is the structure tensor -- it's a symmetric 2x2 matrix for each pixel in the image, based on the covariance of nearby gradients.\n",
    "\n",
    "In other words, it tells us how aligned image gradients are within each region of the image as well as how large they are.\n",
    "\n",
    "Let's dive right in to a farily large standalone example to understand what's going on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio as rio\n",
    "\n",
    "from skimage.feature import structure_tensor, structure_tensor_eigvals\n",
    "\n",
    "from context import data\n",
    "\n",
    "with rio.open(data.naip.lineaments, 'r') as src:\n",
    "    image = src.read()\n",
    "\n",
    "# This assumes a grayscale image. For simplicity, we'll just use RGB mean.\n",
    "data = image.astype(float).mean(axis=0)\n",
    "\n",
    "# Compute the structure tensor. This is basically local gradient similarity.\n",
    "# We're getting three components at each pixel that correspond to a 2x2\n",
    "# symmetric matrix. i.e. [[axx, axy],[axy, ayy]]\n",
    "axx, axy, ayy = structure_tensor(data, sigma=2.5, mode='mirror')\n",
    "\n",
    "# Then we'll compute the eigenvalues of that matrix.\n",
    "v1, v2 = structure_tensor_eigvals(axx, axy, ayy)\n",
    "\n",
    "# And calculate the eigenvector corresponding to the largest eigenvalue.\n",
    "dx, dy = v1 - axx, -axy\n",
    "\n",
    "# We have a vector at each pixel now.  However, we don't really care about all\n",
    "# of them, only those with a large magnitude.  Also, we don't need to worry\n",
    "# about every pixel, as adjacent values are very highly correlated. Therefore,\n",
    "# let's only consider every 10th pixel in each direction.\n",
    "\n",
    "# Top 10th percentile of magnitude\n",
    "mag = np.hypot(dx, dy)\n",
    "selection = mag > np.percentile(mag, 90)\n",
    "\n",
    "# Every 10th pixel (skipping left edge due to boundary effects)\n",
    "ds = np.zeros_like(selection)\n",
    "ds[::10, 10::10] = True\n",
    "selection = ds & selection\n",
    "\n",
    "\n",
    "# Now we'll visualize the selected (large) structure tensor directions both\n",
    "# superimposed on the image and as a rose diagram...\n",
    "fig = plt.figure(constrained_layout=True)\n",
    "ax1 = fig.add_subplot(2, 1, 1)\n",
    "ax2 = fig.add_subplot(2, 1, 2, projection='polar', theta_offset=np.pi/2,\n",
    "                      theta_direction=-1)\n",
    "\n",
    "ax1.imshow(np.moveaxis(image, 0, -1))\n",
    "\n",
    "y, x = np.mgrid[:dx.shape[0], :dx.shape[1]]\n",
    "\n",
    "no_arrow = dict(headwidth=0, headlength=0, headaxislength=0)\n",
    "ax1.quiver(x[selection], y[selection], dx[selection], dy[selection],\n",
    "           angles='xy', units='xy', pivot='middle', color='red', **no_arrow)\n",
    "\n",
    "\n",
    "# We actually want to be perpendictular to the direction of change.. i.e.\n",
    "# we want to point _along_ the lineament. Therefore we'll subtract 90 degrees.\n",
    "# (Could have just gotten the direction of the smaller eigenvector, but we\n",
    "# need to base the magnitude on the largest eigenvector.)\n",
    "angle = np.arctan2(dy[selection], dx[selection]) - np.pi/2\n",
    "ax2.hist(np.concatenate([angle.ravel(), angle.ravel() + np.pi]), bins=120)\n",
    "\n",
    "ax1.set(xticks=[], yticks=[])\n",
    "ax2.set(xticks=[], yticks=[], axisbelow=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thin Section Grain Analysis\n",
    "-----------------------------------------\n",
    "\n",
    "We've done a lot so far with large-scale data.  A key reason to know image processing methods, though, is that they apply to all scales.  Let's shift gears and spend a bit of time working with photomicrographs.  \n",
    "\n",
    "In this case, we'll try to measure the shape preferred orientation of mineral grains in a rock.  We won't go too deeply into the specifics, so let's just try to get an idea of the distribution of grain orientations in our sample.  \n",
    "\n",
    "If you're not familiar with optical petrology, we're looking at a slide with a very thin (~30 microns) slice of rock attached to it.  It's common to use both plane polarized light and cross polarized light.  You can compare the change in appearance in the figure below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio as rio\n",
    "import skimage.io\n",
    "\n",
    "from context import data\n",
    "from context import utils\n",
    "\n",
    "xpl_rgb = skimage.io.imread(data.bgs_rock.amphibolite_xpl)\n",
    "ppl_rgb = skimage.io.imread(data.bgs_rock.amphibolite_ppl)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "im_ppl = ax.imshow(ppl_rgb, label='Plane Polarized')\n",
    "im_xpl = ax.imshow(xpl_rgb, label='Cross Polarized')\n",
    "ax.set(xticks=[], yticks=[])\n",
    "utils.Toggler(im_xpl).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the (poor) 40,000 foot (yes, I know... Units...) overview: \n",
    "\n",
    "Minerals with a high degree of birefringence (controlled by anisotropy in the speed of light through the crystal) can appear in a range of \"gaudy\" colors under cross polars.  The color depends on the orientation of the crystal, so an individual color doesn't tell us much. The range of colors for a certain mineral is a useful indicator, though.  \n",
    "\n",
    "In this case, we don't really care what the colors mean, we only want to use the color to distinguish one grain from another adjacent grain.  We're interested in being able to identify distinct grains because we're interested in the shape of each grain.\n",
    "\n",
    "If you look closely, you'll see that the colors within a grain aren't constant. They vary a bit. There are also lots of small mineral grains and inclusions within larger grains that we don't need to worry as much about for this analysis. Ideally we'd try to separate them, but it's okay if we don't. We mostly want to look at the orientation of the largest grains.\n",
    "\n",
    "Okay. So we want to identify distinct grains.  In image processing terms, we want to segment the image. We'll do this by finding regions with similar colors.  These \"regions of a similar color\" often are referred to as \"superpixels\" in image processing terms.  There are a _huge_ variety of methods to do this. Some use gradients to define \"watersheds\" of similar color, some use clustering, and you can even apply more flexible methods like a trained CNN to do this.\n",
    "\n",
    "Let's apply a fairly well-known and widely used \"superpixel\" method: Simple Linear Iterative Clustering (SLIC).  It's based on [K-means clustering](https://en.wikipedia.org/wiki/K-means_clustering), which is a common method to find groups of similar data.  In this case,  it's clustering spatially as well as in color-space.  It's actually using X,Y coordinates directly in the clustering and working in 5 dimensions (RGB + XY).\n",
    "\n",
    "To make things a bit easier on ourselves, let's just work with the center of the image (avoids needing to separate out the background that's not part of the thin section):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "xpl_rgb = xpl_rgb[500:3000, 1000:4000, :]\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(xpl_rgb)\n",
    "ax.set(xticks=[], yticks=[])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's segment the image using SLIC and play around with the parameters a bit..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage.segmentation\n",
    "\n",
    "grains = skimage.segmentation.slic(xpl_rgb, sigma=0.5, multichannel=True,\n",
    "                                   n_segments=1500, compactness=0.1)\n",
    "\n",
    "# It's hard to color each grain with a unique color, so we'll show boundaries\n",
    "# in yellow instead of coloring them like we did before.\n",
    "overlay = skimage.segmentation.mark_boundaries(xpl_rgb, grains, start_label=1)\n",
    "\n",
    "fig, ax = plt.subplots(constrained_layout=True)\n",
    "ax.imshow(overlay)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a more compelete example where we extract properties (such as long axis and short axis) of each segmented region to look at the distribution of grain orientations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import skimage.io\n",
    "import skimage.segmentation\n",
    "import skimage.measure\n",
    "import scipy.ndimage\n",
    "\n",
    "from context import data\n",
    "from context import utils\n",
    "\n",
    "# Amphibolite under cross polars (from BGS, see data/bgs_rock/README.md)\n",
    "xpl_rgb = skimage.io.imread(data.bgs_rock.amphibolite_xpl)\n",
    "\n",
    "# Let's use the center of the image to avoid needing to worry about the edges.\n",
    "xpl_rgb = xpl_rgb[500:3000, 1000:4000, :]\n",
    "\n",
    "# This attempts to group locally similar colors. It's kmeans in 5 dimensions\n",
    "# (RGB + XY).  N_segments and compactness are the main \"knobs\" to turn.\n",
    "grains = skimage.segmentation.slic(xpl_rgb, sigma=0.5, multichannel=True,\n",
    "                                   n_segments=1500, compactness=0.1)\n",
    "\n",
    "# It's hard to color each grain with a unique color, so we'll show boundaries\n",
    "# in yellow instead of coloring them like we did before.\n",
    "overlay = skimage.segmentation.mark_boundaries(xpl_rgb, grains)\n",
    "\n",
    "# Now let's extract information about each individual grain we've classified.\n",
    "# In this case, we're only interested in orientation, but there's a lot more\n",
    "# we could extract.\n",
    "info = skimage.measure.regionprops(grains)\n",
    "\n",
    "# And calculate the orientation of the long axis of each grain...\n",
    "angles = []\n",
    "for item in info:\n",
    "    cov = item['inertia_tensor']\n",
    "    azi = np.degrees(np.arctan2((-2 * cov[0, 1]), (cov[0,0] - cov[1,1])))\n",
    "    angles.append(azi)\n",
    "\n",
    "# Make bidirectional (quick hack for plotting)\n",
    "angles = angles + [x + 180 for x in angles]\n",
    "\n",
    "# Now display the segmentation and a rose diagram\n",
    "fig = plt.figure(constrained_layout=True)\n",
    "ax1 = fig.add_subplot(1, 2, 1)\n",
    "ax2 = fig.add_subplot(1, 2, 2, projection='polar', theta_offset=np.pi/2,\n",
    "                      theta_direction=-1)\n",
    "ax1.imshow(overlay)\n",
    "ax2.hist(np.radians(angles), bins=60)\n",
    "\n",
    "ax1.set(xticks=[], yticks=[])\n",
    "ax2.set(xticklabels=[], yticklabels=[], axisbelow=True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wrapping it all up\n",
    "\n",
    "That's all for now, folks!  Hopefully you can see ways to apply some of this to problems you're actively working on.  There are a ton of very powerful methods exposed in common python image processing libraries, and it's easy to get a bit lost in the huge variety of options.  Hopefully this has given you enough of an understanding of common methods to start exploring on your own.  There's a lot out there there's very "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}