Alignment-RigRelatives.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Passive Image Alignment\n",
    "\n",
    "This workbook shows a completely passive alignment function using only the RigRelatives present in image metadata.  Older versions of firmware may not include these tags. If your images don't, check the MicaSense [support site](https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields) for some tips on how to update your camera firmware to have them, as well as how to [add them to datasets](https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support) prior to the update. \n",
    "\n",
    "While this method doesn't provide perfect alignment, it can be fast and very useful for visualization of images when processing power is limited or speed is more important than alignment quality.\n",
    "\n",
    "## Opening Images\n",
    "\n",
    "As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize the 5 bands of a RedEdge capture.\n",
    "\n",
    "First, we'll load the `autoreload` extension.  This lets us change underlying code (such as library functions) without having to reload the entire workbook and kernel. This is useful in this workbook because the cell that runs the alignment can take a long time to run, so with `autoreload` extension we can change external code for analysis and visualization without needing to re-compute the alignments each time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import os, glob\n",
    "import micasense.capture as capture\n",
    "%matplotlib inline\n",
    "from pathlib import Path\n",
    "\n",
    "panelNames = None\n",
    "\n",
    "# This is an Altum image with RigRelatives and a thermal band\n",
    "\n",
    "imagePath = Path(\"./data/ALTUM\")\n",
    "\n",
    "# these will return lists of image paths as strings \n",
    "imageNames = list(imagePath.glob('IMG_0021_*.tif'))\n",
    "imageNames = [x.as_posix() for x in imageNames]\n",
    "\n",
    "panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n",
    "panelNames = [x.as_posix() for x in panelNames]\n",
    "\n",
    "if panelNames is not None:\n",
    "    panelCap = capture.Capture.from_filelist(panelNames)\n",
    "else:\n",
    "    panelCap = None\n",
    "\n",
    "capture = capture.Capture.from_filelist(imageNames)\n",
    "\n",
    "for img in capture.images:\n",
    "    if img.rig_relatives is None:\n",
    "        raise ValueError(\"Images must have RigRelatives tags set which requires updated firmware and calibration. See the links in text above\")\n",
    "\n",
    "if panelCap is not None:\n",
    "    if panelCap.panel_albedo() is not None:\n",
    "        panel_reflectance_by_band = panelCap.panel_albedo()\n",
    "    else:\n",
    "        panel_reflectance_by_band = [0.49, 0.49, 0.49, 0.49, 0.49] #RedEdge band_index order\n",
    "    panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)    \n",
    "    img_type = \"reflectance\"\n",
    "    capture.plot_undistorted_reflectance(panel_irradiance)\n",
    "else:\n",
    "    if False: #capture.dls_present():\n",
    "        img_type='reflectance'\n",
    "        capture.plot_undistorted_reflectance(capture.dls_irradiance())\n",
    "    else:\n",
    "        img_type = \"radiance\"\n",
    "        capture.plot_undistorted_radiance()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alignment - Rig Relatives\n",
    "\n",
    "For images with RigRelative tags present, we can find a rough alignment using only the built in relatives. These can be good for quick visualizatoins.  For better results, use an algorithm like that in the other image alignment tutorial to sweeten the alignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import micasense.imageutils as imageutils\n",
    "import micasense.plotutils as plotutils\n",
    "\n",
    "reference_band = 5\n",
    "warp_mode = cv2.MOTION_HOMOGRAPHY\n",
    "warp_matrices = capture.get_warp_matrices(ref_index=reference_band)\n",
    "\n",
    "cropped_dimensions,edges = imageutils.find_crop_bounds(capture,warp_matrices,reference_band=reference_band)\n",
    "im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, reference_band, img_type=img_type)\n",
    "\n",
    "print(im_aligned.shape)\n",
    "\n",
    "print(\"warp_matrices={}\".format(warp_matrices))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Aligned Images\n",
    "\n",
    "Once the transformation has been found, it can be verified by composting the aligned images to check alignment. The image 'stack' containing all bands can also be exported to a multi-band TIFF file for viewing in extrernal software such as QGIS.  Usef ul componsites are a naturally colored RGB as well as color infrared, or CIR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# figsize=(30,23) # use this size for full-image-resolution display\n",
    "figsize=(16,13)   # use this size for export-sized display\n",
    "\n",
    "rgb_band_indices = [2,1,0]\n",
    "cir_band_indices = [3,2,1]\n",
    "\n",
    "# Create an empty normalized stack for viewing\n",
    "im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32 )\n",
    "\n",
    "im_min = np.percentile(im_aligned[:,:,0:2].flatten(),  0.1)  # modify with these percentilse to adjust contrast\n",
    "im_max = np.percentile(im_aligned[:,:,0:2].flatten(), 99.9)  # for many images, 0.5 and 99.5 are good values\n",
    "\n",
    "for i in range(0,im_aligned.shape[2]):\n",
    "    if img_type == 'reflectance':\n",
    "        # for reflectance images we maintain white-balance by applying the same display scaling to all bands\n",
    "        im_display[:,:,i] =  imageutils.normalize(im_aligned[:,:,i], im_min, im_max)\n",
    "    elif img_type == 'radiance':\n",
    "        # for radiance images we do an auto white balance since we don't know the input light spectrum by\n",
    "        # stretching each display band histogram to it's own min and max\n",
    "        im_display[:,:,i] =  imageutils.normalize(im_aligned[:,:,i])\n",
    "\n",
    "rgb = im_display[:,:,rgb_band_indices]\n",
    "# for cir false color imagery, we normalize the NIR,R,G bands within themselves, which provides\n",
    "# the classical CIR rendering where plants are red and soil takes on a blue tint\n",
    "for i in cir_band_indices:\n",
    "    im_display[:,:,i] =  imageutils.normalize(im_aligned[:,:,i])\n",
    "\n",
    "cir = im_display[:,:,cir_band_indices]\n",
    "fig, axes = plt.subplots(1, 2, figsize=figsize)\n",
    "axes[0].set_title(\"Red-Green-Blue Composite\")\n",
    "axes[0].imshow(rgb)\n",
    "axes[1].set_title(\"Color Infrared (CIR) Composite\")\n",
    "axes[1].imshow(cir)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Enhancement\n",
    "\n",
    "There are many techniques for image enhancement, but one which is commonly used to improve the visual sharpness of imagery is the unsharp mask.  Here we apply an unsharp mask to the RGB image to improve the visualization, and then apply a gamma curve to make the darkest areas brighter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create an enhanced version of the RGB render using an unsharp mask\n",
    "gaussian_rgb = cv2.GaussianBlur(rgb, (9,9), 10.0)\n",
    "gaussian_rgb[gaussian_rgb<0] = 0\n",
    "gaussian_rgb[gaussian_rgb>1] = 1\n",
    "unsharp_rgb = cv2.addWeighted(rgb, 1.5, gaussian_rgb, -0.5, 0)\n",
    "unsharp_rgb[unsharp_rgb<0] = 0\n",
    "unsharp_rgb[unsharp_rgb>1] = 1\n",
    "\n",
    "# Apply a gamma correction to make the render appear closer to what our eyes would see\n",
    "gamma = 1.4\n",
    "gamma_corr_rgb = unsharp_rgb**(1.0/gamma)\n",
    "fig = plt.figure(figsize=figsize)\n",
    "plt.imshow(gamma_corr_rgb, aspect='equal')\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rendered Image output\n",
    "\n",
    "We can output the image to a PNG or JPEG file for viewing. This can also be useful in creating thumbnails of captures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio\n",
    "imtype = 'png' # or 'jpg'\n",
    "imageio.imwrite('rgb.'+imtype, (255*gamma_corr_rgb).astype('uint8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stack Export\n",
    "\n",
    "We can export the image easily stacks using the `gdal` library (http://www.glal.org). Once exported, these image stacks can be opened in software such as QGIS and raster operations such as NDVI or NDRE computation can be done in that software.  At this time the stacks don't include any geographic information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal, gdal_array\n",
    "rows, cols, bands = im_display.shape\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "filename = \"bgrne\" #blue,green,red,nir,redEdge\n",
    "\n",
    "if im_aligned.shape[2] == 6:\n",
    "    filename = filename + \"t\" #thermal\n",
    "outRaster = driver.Create(filename+\".tiff\", cols, rows, im_aligned.shape[2], gdal.GDT_UInt16)\n",
    "\n",
    "normalize = (img_type == 'radiance') # normalize radiance images to fit with in UInt16\n",
    "\n",
    "# Output a 'stack' in the same band order as RedEdge/Alutm\n",
    "# Blue,Green,Red,NIR,RedEdge[,Thermal]\n",
    "# reflectance stacks are output with 32768=100% reflectance to provide some overhead for specular reflections\n",
    "# radiance stacks are output with 65535=100% radiance to provide some overhead for specular reflections\n",
    "\n",
    "# NOTE: NIR and RedEdge are not in wavelength order!\n",
    "\n",
    "multispec_min = np.min(im_aligned[:,:,1:5])\n",
    "multispec_max = np.max(im_aligned[:,:,1:5])\n",
    "\n",
    "for i in range(0,5):\n",
    "    outband = outRaster.GetRasterBand(i+1)\n",
    "    if normalize:\n",
    "        outdata = imageutils.normalize(im_aligned[:,:,i],multispec_min,multispec_max)\n",
    "    else:\n",
    "        outdata = im_aligned[:,:,i]\n",
    "        outdata[outdata<0] = 0\n",
    "        outdata[outdata>2] = 2\n",
    "    \n",
    "    outdata = outdata*32767\n",
    "    outdata[outdata<0] = 0\n",
    "    outdata[outdata>65535] = 65535\n",
    "    outband.WriteArray(outdata)\n",
    "    outband.FlushCache()\n",
    "\n",
    "if im_aligned.shape[2] == 6:\n",
    "    outband = outRaster.GetRasterBand(6)\n",
    "    outdata = im_aligned[:,:,5] * 100 # scale to centi-C to fit into uint16\n",
    "    outdata[outdata<0] = 0\n",
    "    outdata[outdata>65535] = 65535\n",
    "    outband.WriteArray(outdata)\n",
    "    outband.FlushCache()\n",
    "outRaster = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes on Alignment and Stack Usage\n",
    "\n",
    "\"Stacks\" as described above are useful in a number of processing cases.  For example, at the time of this writing, many photogrammetry suites could import and process stack files without significantly impacting the radiometric processing which has already been accomplished.  \n",
    "\n",
    "Running photogrammetry on stack files instead of raw image files has both advantages and drawbacks. The primary advantage has been found to be an increase in processing speed and a reduction in program memory usage. As the photogrammetric workflow generally operates on luminance images and may not use color information, stacked images may require similar resources and be processed at a similar speed as single-band images.  This is because one band of the stack can be used to generate the matching feature space while the others are ignored for matching purposes. This reduces the feature space 5-fold over matching using all images separately.\n",
    "\n",
    "One disadvantage is that stacking images outside of the photogrammetric workflow may result in poor image matching.  The RedEdge is known to have stable lens characteristics over the course of normal operation, but variations in temperature or impacts to the camera through handling or rough landings may change the image alignment parameters.  For this reason, we recommend finding a matching transformation for each flight (each take-off and landing).  Alignment transformations from multiple images within a flight can be compared to find the best transformation to apply to the set of the flight.  While not described or supported in this generic implementation, some matching algorithms can use a \"seed\" value as a starting point to speed up matching.  For most cases, this seed could be the transformation found in a previous flight, or another source of a known good transformation.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NDVI Computation\n",
    "\n",
    "For raw index computation on single images, the numpy package provides a simple way to do math and simple visualizatoin on images.  Below, we compute and visualize an image histogram and then use that to pick a colormap range for visualizing the NDVI of an image. \n",
    "\n",
    "### Plant Classification\n",
    "\n",
    "After computing the NDVI and prior to displaying it, we use a very rudimentary method for focusing on the plants and removing the soil and shadow information from our images and histograms. Below we remove non-plant pixels by setting to zero any pixels in the image where the NIR reflectance is less than 20%.  This helps to ensure that the NDVI and NDRE histograms aren't skewed substantially by soil noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from micasense import plotutils\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation\n",
    "\n",
    "# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands\n",
    "ndvi = (im_aligned[:,:,3] - im_aligned[:,:,2]) / (im_aligned[:,:,3] + im_aligned[:,:,2])\n",
    "\n",
    "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n",
    "if img_type == 'reflectance':\n",
    "    ndvi = np.ma.masked_where(im_aligned[:,:,3] < 0.20, ndvi) \n",
    "elif img_type == 'radiance':\n",
    "    lower_pct_radiance = np.percentile(im_aligned[:,:,3],  10.0)\n",
    "    ndvi = np.ma.masked_where(im_aligned[:,:,3] < lower_pct_radiance, ndvi) \n",
    "    \n",
    "# Compute and display a histogram\n",
    "ndvi_hist_min = np.min(ndvi)\n",
    "ndvi_hist_max = np.max(ndvi)\n",
    "fig, axis = plt.subplots(1, 1, figsize=(10,4))\n",
    "axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))\n",
    "plt.title(\"NDVI Histogram\")\n",
    "plt.show()\n",
    "\n",
    "min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values\n",
    "#min_display_ndvi = np.percentile(ndvi.flatten(),  5.0)  # modify with these percentilse to adjust contrast\n",
    "max_display_ndvi = np.percentile(ndvi.flatten(), 99.5)  # for many images, 0.5 and 99.5 are good values\n",
    "masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)\n",
    "\n",
    "#reduce the figure size to account for colorbar\n",
    "figsize=np.asarray(figsize) - np.array([3,2])\n",
    "\n",
    "#plot NDVI over an RGB basemap, with a colorbar showing the NDVI scale\n",
    "fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb, \n",
    "                                    masked_ndvi, \n",
    "                                    figsize = figsize, \n",
    "                                    title = 'NDVI filtered to only plants over RGB base layer',\n",
    "                                    vmin = min_display_ndvi,\n",
    "                                    vmax = max_display_ndvi)\n",
    "fig.savefig('ndvi_over_rgb.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NDRE Computation\n",
    "\n",
    "In the same manner, we can compute, filter, and display another index useful for the RedEdge camera, the Normalized Difference Red Edge (NDRE) index.  We also filter out shadows and soil to ensure our display focuses only on the plant health."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n",
    "ndre = (im_aligned[:,:,3] - im_aligned[:,:,4]) / (im_aligned[:,:,3] + im_aligned[:,:,4])\n",
    "\n",
    "# Mask areas with shadows and low NDVI to remove soil\n",
    "masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)\n",
    "\n",
    "# Compute a histogram\n",
    "ndre_hist_min = np.min(masked_ndre)\n",
    "ndre_hist_max = np.max(masked_ndre)\n",
    "fig, axis = plt.subplots(1, 1, figsize=(10,4))\n",
    "axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))\n",
    "plt.title(\"NDRE Histogram (filtered to only plants)\")\n",
    "plt.show()\n",
    "\n",
    "min_display_ndre = np.percentile(masked_ndre, 5)\n",
    "max_display_ndre = np.percentile(masked_ndre, 99.5)\n",
    "\n",
    "fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb, \n",
    "                                    masked_ndre, \n",
    "                                    figsize=figsize, \n",
    "                                    title='NDRE filtered to only plants over RGB base layer',\n",
    "                                    vmin=min_display_ndre,vmax=max_display_ndre)\n",
    "fig.savefig('ndre_over_rgb.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thermal Imagery\n",
    "\n",
    "If our image is from an Altum and includes a thermal band, we can display the re-sampled and aligned thermal data over the RGB data to maintain the context of the thermal information.\n",
    "\n",
    "In the image below, it's very clear based on the average temperature where the soil is wet and dry, and even in the middle of the road we can find some wetter areas of soil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if im_aligned.shape[2] >= 5:\n",
    "\n",
    "    # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n",
    "    masked_thermal = im_aligned[:,:,5]\n",
    "    # Alternatively we can mask the thermal only to plants here, which is useful for large contiguous areas\n",
    "    # masked_thermal = np.ma.masked_where(ndvi < 0.45, im_aligned[:,:,5])\n",
    "\n",
    "\n",
    "    # Compute a histogram\n",
    "    fig, axis = plt.subplots(1, 1, figsize=(10,4))\n",
    "    axis.hist(masked_thermal.ravel(), bins=512, range=(np.min(masked_thermal), np.max(masked_thermal)))\n",
    "    plt.title(\"Thermal Histogram\")\n",
    "    plt.show()\n",
    "\n",
    "    min_display_therm = np.percentile(masked_thermal, 1)\n",
    "    max_display_therm = np.percentile(masked_thermal, 99)\n",
    "\n",
    "    fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,\n",
    "                                        masked_thermal, \n",
    "                                        figsize=figsize, \n",
    "                                        title='Temperature over True Color',\n",
    "                                        vmin=min_display_therm,vmax=max_display_therm,\n",
    "                                        overlay_alpha=0.25,\n",
    "                                        overlay_colormap='jet',\n",
    "                                        overlay_steps=16,\n",
    "                                        display_contours=True,\n",
    "                                        contour_steps=16,\n",
    "                                        contour_alpha=.4,\n",
    "                                        contour_fmt=\"%.0fC\")\n",
    "    fig.savefig('thermal_over_rgb.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Copyright (c) 2017-2019 MicaSense, Inc.  For licensing information see the [project git repository](https://github.com/micasense/imageprocessing)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.12"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}