diff --git a/Alignment-10Band.ipynb b/Alignment-10Band.ipynb deleted file mode 100644 index 87d05a87..00000000 --- a/Alignment-10Band.ipynb +++ /dev/null @@ -1,542 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Active Image Alignment - 10-band\n", - "\n", - "Note: This notebook specifically demonstrates alignment and stack creation for 10-band data. For alignment of 5-band (RedEdge) and 6-band (Altum) data, pr more details about alignment in general, see the \"Active Image Alignment\" notebook.\n", - "\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 10 bands of a RedEdge-MX Dual (Red+Blue) 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 update the code after the alignment step 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": {}, - "outputs": [], - "source": [ - "import os, glob\n", - "import micasense.capture as capture\n", - "%matplotlib inline\n", - "\n", - "panelNames = None\n", - "paneCap = None\n", - "\n", - "imagePath = os.path.join('.','data','10BANDSET','000')\n", - "imageNames = glob.glob(os.path.join(imagePath,'IMG_0431_*.tif'))\n", - "panelNames = glob.glob(os.path.join(imagePath,'IMG_0000_*.tif'))\n", - "\n", - "# Allow this code to align both radiance and reflectance images; bu excluding\n", - "# a definition for panelNames above, radiance images will be used\n", - "# For panel images, efforts will be made to automatically extract the panel information\n", - "# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance\n", - "# will need to be set in the panel_reflectance_by_band variable.\n", - "# Note: radiance images will not be used to properly create NDVI/NDRE images below.\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", - "if panelCap is not None:\n", - " if panelCap.panel_albedo() is not None:\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " raise IOError(\"Comment this lne and set panel_reflectance_by_band here\")\n", - " panel_reflectance_by_band = [0.65]*len(imageNames)\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 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": [ - "## Unwarp and Align\n", - "\n", - "Alignment is a three step process:\n", - "\n", - "1. Images are unwarped using the built-in lens calibration\n", - "1. A transformation is found to align each band to a common band\n", - "1. The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", - "\n", - "We provide utilities to find the alignement transformations within a single capture. Our experience shows that once a good alignmnet transformation is found, it tends to be stable over a flight and, in most cases, over many flights. The transformation may change if the camera undergoes a shock event (such as a hard landing or drop) or if the temperature changes substantially between flights. In these events a new transformation may need to be found.\n", - "\n", - "Further, since this approach finds a 2-dimensional (affine) transformation between images, it won't work when the parallax between bands results in a 3-dimensional depth field. This can happen if very close to the target or when targets are visible at significantly different ranges, such as a nearby tree or building against a background much farther way. In these cases it will be necessary to use photogrammetry techniques to find a 3-dimensional mapping between images.\n", - "\n", - "For best alignment results it's good to select a capture which has features which visible in all bands. Man-made objects such as cars, roads, and buildings tend to work very well, while captures of only repeating crop rows tend to work poorly. Remember, once a good transformation has been found for flight, it can be generally be applied across all of the images.\n", - "\n", - "It's also good to use an image for alignment which is taken near the same level above ground as the rest of the flights. Above approximately 35m AGL, the alignement will be consistent. However, if images taken closer to the ground are used, such as panel images, the same alignment transformation will not work for the flight data. " - ] - }, - { - "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", - "## Alignment settings\n", - "match_index = 4 # Index of the band, here we use green\n", - "max_alignment_iterations = 20\n", - "warp_mode = cv2.MOTION_HOMOGRAPHY # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", - "pyramid_levels = 3 # for 10-band imagery we use a 3-level pyramid. In some cases\n", - "\n", - "print(\"Alinging images. Depending on settings this can take from a few seconds to many minutes\")\n", - "# Can potentially increase max_iterations for better results, but longer runtimes\n", - "warp_matrices, alignment_pairs = imageutils.align_capture(capture,\n", - " ref_index = match_index,\n", - " max_iterations = max_alignment_iterations,\n", - " warp_mode = warp_mode,\n", - " pyramid_levels = pyramid_levels)\n", - "\n", - "print(\"Finished Aligning, warp matrices={}\".format(warp_matrices))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Crop Aligned Images\n", - "After finding image alignments we may need to remove pixels around the edges which aren't present in every image in the capture. To do this we use the affine transforms found above and the image distortions from the image metadata. OpenCV provides a couple of handy helpers for this task in the `cv2.undistortPoints()` and `cv2.transform()` methods. These methods takes a set of pixel coordinates and apply our undistortion matrix and our affine transform, respectively. So, just as we did when registering the images, we first apply the undistortion process the coordinates of the image borders, then we apply the affine transformation to that result. The resulting pixel coordinates tell us where the image borders end up after this pair of transformations, and we can then crop the resultant image to these coordinates." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cropped_dimensions, edges = imageutils.find_crop_bounds(capture, warp_matrices, warp_mode=warp_mode)\n", - "im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, match_index, img_type=img_type)" - ] - }, - { - "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": { - "scrolled": false - }, - "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 = [capture.band_names().index('Red'),capture.band_names().index('Green'),capture.band_names().index('Blue-444')]\n", - "cir_band_indices = [capture.band_names().index('NIR'),capture.band_names().index('Red'),capture.band_names().index('Green')]\n", - "\n", - "# Create a 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[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast\n", - "im_max = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values\n", - "\n", - "# for rgb true color, we use the same min and max scaling across the 3 bands to \n", - "# maintain the \"white balance\" of the calibrated image\n", - "for i in rgb_band_indices:\n", - " im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)\n", - "\n", - "rgb = im_display[:,:,rgb_band_indices]\n", - "\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": {}, - "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": [ - "## Image Export\n", - "\n", - "Composite images can be exported to JPEG or PNG format using the `imageio` package. These images may be useful for visualization or thumbnailing, and creating RGB thumbnails of a set of images can provide a convenient way to browse the imagery in a more visually appealing way that browsing the raw imagery. " - ] - }, - { - "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'))\n", - "imageio.imwrite('cir.'+imtype, (255*cir).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 = \"bbggrreeen\" #blue,blue,green,green,red,red,rededge,rededge,rededge,nir\n", - "\n", - "sort_by_wavelength = True # set to false if you want stacks in camera-band-index order\n", - " \n", - "outRaster = driver.Create(filename+\".tiff\", cols, rows, bands, gdal.GDT_UInt16, options = [ 'INTERLEAVE=BAND','COMPRESS=DEFLATE' ])\n", - "try:\n", - " if outRaster is None:\n", - " raise IOError(\"could not load gdal GeoTiff driver\")\n", - "\n", - " if sort_by_wavelength:\n", - " eo_list = list(np.argsort(capture.center_wavelengths()))\n", - " else:\n", - " eo_list = capture.eo_indices()\n", - "\n", - " for outband,inband in enumerate(eo_list):\n", - " outband = outRaster.GetRasterBand(outband+1)\n", - " outdata = im_aligned[:,:,inband]\n", - " outdata[outdata<0] = 0\n", - " outdata[outdata>2] = 2 #limit reflectance data to 200% to allow some specular reflections\n", - " outband.WriteArray(outdata*32768) # scale reflectance images so 100% = 32768\n", - " outband.FlushCache()\n", - "\n", - " for outband,inband in enumerate(capture.lw_indices()):\n", - " outband = outRaster.GetRasterBand(len(eo_list)+outband+1)\n", - " outdata = (im_aligned[:,:,inband]+273.15) * 100 # scale data from float degC to back to centi-Kelvin to fit into uint16\n", - " outdata[outdata<0] = 0\n", - " outdata[outdata>65535] = 65535\n", - " outband.WriteArray(outdata)\n", - " outband.FlushCache()\n", - "finally:\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", - "nir_band = [name.lower() for name in capture.band_names()].index('nir')\n", - "red_band = [name.lower() for name in capture.band_names()].index('red')\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[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])\n", - "\n", - "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", - "if img_type == 'reflectance':\n", - " ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < 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[:,:,nir_band] < 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", - "rededge_band = [name.lower() for name in capture.band_names()].index('red edge')\n", - "ndre = (im_aligned[:,:,nir_band] - im_aligned[:,:,rededge_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,rededge_band])\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": [ - "## Red vs NIR Reflectance\n", - "\n", - "Finally, we show a classic agricultural remote sensing output in the tassled cap plot. This plot can be useful for visualizing row crops and plots the Red Reflectance channel on the X-axis against the NIR reflectance channel on the Y-axis. This plot also clearly shows the line of the soil in that space. The tassled cap view isn't very useful for this arid data set; however, we can see the \"badge of trees\" of high NIR reflectance and relatively low red reflectance. This provides an example of one of the uses of aligned images for single capture analysis." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "x_band = red_band\n", - "y_band = nir_band\n", - "x_max = np.max(im_aligned[:,:,x_band])\n", - "y_max = np.max(im_aligned[:,:,y_band])\n", - "\n", - "fig = plt.figure(figsize=(12,12))\n", - "plt.hexbin(im_aligned[:,:,x_band],im_aligned[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))\n", - "ax = fig.gca()\n", - "ax.set_xlim([0,x_max])\n", - "ax.set_ylim([0,y_max])\n", - "plt.xlabel(\"{} Reflectance\".format(capture.band_names()[x_band]))\n", - "plt.ylabel(\"{} Reflectance\".format(capture.band_names()[y_band]))\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Print warp_matrices for usage elsewhere, such as Batch Processing\n", - "Last, we output the warp_matrices that we got for this image stack for usage elsewhere. Currently these can be used in the `Batch Processing.ipynb` notebook to save reflectance-compensated stacks of images to a directory." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "print(warp_matrices)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "---\n", - "Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the [project git repository](https://github.com/micasense/imageprocessing)" - ] - } - ], - "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.7.3" - }, - "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 -} diff --git a/Alignment-RigRelatives.ipynb b/Alignment-RigRelatives.ipynb deleted file mode 100644 index 215318b3..00000000 --- a/Alignment-RigRelatives.ipynb +++ /dev/null @@ -1,512 +0,0 @@ -{ - "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", - "\n", - "panelNames = None\n", - "\n", - "# This is an altum image with RigRelatives and a thermal band\n", - "imagePath = os.path.join('.','data','ALTUM1SET','000')\n", - "imageNames = glob.glob(os.path.join(imagePath,'IMG_0008_*.tif'))\n", - "panelNames = glob.glob(os.path.join(imagePath,'IMG_0000_*.tif'))\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.67, 0.69, 0.68, 0.61, 0.67] #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", - "warp_mode = cv2.MOTION_HOMOGRAPHY\n", - "warp_matrices = capture.get_warp_matrices()\n", - "\n", - "cropped_dimensions,edges = imageutils.find_crop_bounds(capture,warp_matrices)\n", - "im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, None, img_type=img_type)\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],capture.num_bands+1), 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", - "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.3" - }, - "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 -} diff --git a/Alignment.ipynb b/Alignment.ipynb deleted file mode 100644 index 9e416d38..00000000 --- a/Alignment.ipynb +++ /dev/null @@ -1,620 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Active Image Alignment\n", - "\n", - "In most use cases, each band of a multispectral capture must be aligned with the other bands in order to create meaningful data. In this tutorial, we show how to align the band to each other using open source OpenCV utilities. \n", - "\n", - "Image alignment allows the combination of images into true-color (RGB) and false color (such as CIR) composites, useful for scouting using single images as well as for display and management uses. In addition to composite images, alignment allows the calculation of pixel-accurate indices such as NDVI or NDRE at the single image level which can be very useful for applications like plant counting and coverage estimations, where mosaicing artifacts may otherwise skew analysis results. \n", - "\n", - "The image alignment method described below tends to work well on images with abundant image features, or areas of significant contrast. Cars, buildings, parking lots, and roads tend to provide the best results. This approach may not work well on images which contain few features or very repetitive features, such as full canopy row crops or fields of repetitive small crops such lettuce or strawberries. We will disscuss more about the advantages and disadvantages of these methods below.\n", - "\n", - "The functions behind this alignment process can work with most versions of RedEdge and Altum firmware. They will work best with versions above 3.2.0 which include the \"RigRelatives\" tags. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.\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 update the code after the alignment step 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": {}, - "outputs": [], - "source": [ - "import os, glob\n", - "import micasense.capture as capture\n", - "%matplotlib inline\n", - "\n", - "panelNames = None\n", - "\n", - "# # This is an older RedEdge image without RigRelatives\n", - "# imagePath = os.path.join(os.path.abspath('.'),'data','0000SET','000')\n", - "# imageNames = glob.glob(os.path.join(imagePath,'IMG_0001_*.tif'))\n", - "# panelNames = glob.glob(os.path.join(imagePath,'IMG_0000_*.tif'))\n", - "\n", - "# # Image from the example RedEdge imageSet (see the ImageSet notebook) without RigRelatives.\n", - "# imagePath = os.path.expanduser(os.path.join('~','Downloads','RedEdgeImageSet','0000SET'))\n", - "# imageNames = glob.glob(os.path.join(imagePath,'000','IMG_0013_*.tif'))\n", - "# panelNames = glob.glob(os.path.join(imagePath,'000','IMG_0000_*.tif'))\n", - "\n", - "# This is an altum image with RigRelatives and a thermal band\n", - "imagePath = os.path.join('.','data','ALTUM1SET','000')\n", - "imageNames = glob.glob(os.path.join(imagePath,'IMG_0245_*.tif'))\n", - "panelNames = glob.glob(os.path.join(imagePath,'IMG_0000_*.tif'))\n", - "\n", - "\n", - "# Allow this code to align both radiance and reflectance images; bu excluding\n", - "# a definition for panelNames above, radiance images will be used\n", - "# For panel images, efforts will be made to automatically extract the panel information\n", - "# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance\n", - "# will need to be set in the panel_reflectance_by_band variable.\n", - "# Note: radiance images will not be used to properly create NDVI/NDRE images below.\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", - "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.67, 0.69, 0.68, 0.61, 0.67] #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 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": [ - "## Unwarp and Align\n", - "\n", - "Alignment is a three step process:\n", - "\n", - "1. Images are unwarped using the built-in lens calibration\n", - "1. A transformation is found to align each band to a common band\n", - "1. The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", - "\n", - "We provide utilities to find the alignement transformations within a single capture. Our experience shows that once a good alignmnet transformation is found, it tends to be stable over a flight and, in most cases, over many flights. The transformation may change if the camera undergoes a shock event (such as a hard landing or drop) or if the temperature changes substantially between flights. In these events a new transformation may need to be found.\n", - "\n", - "Further, since this approach finds a 2-dimensional (affine) transformation between images, it won't work when the parallax between bands results in a 3-dimensional depth field. This can happen if very close to the target or when targets are visible at significantly different ranges, such as a nearby tree or building against a background much farther way. In these cases it will be necessary to use photogrammetry techniques to find a 3-dimensional mapping between images.\n", - "\n", - "For best alignment results it's good to select a capture which has features which visible in all bands. Man-made objects such as cars, roads, and buildings tend to work very well, while captures of only repeating crop rows tend to work poorly. Remember, once a good transformation has been found for flight, it can be generally be applied across all of the images.\n", - "\n", - "It's also good to use an image for alignment which is taken near the same level above ground as the rest of the flights. Above approximately 35m AGL, the alignement will be consistent. However, if images taken closer to the ground are used, such as panel images, the same alignment transformation will not work for the flight data. " - ] - }, - { - "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", - "## Alignment settings\n", - "match_index = 1 # Index of the band \n", - "max_alignment_iterations = 10\n", - "warp_mode = cv2.MOTION_HOMOGRAPHY # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", - "pyramid_levels = 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment\n", - "\n", - "print(\"Alinging images. Depending on settings this can take from a few seconds to many minutes\")\n", - "# Can potentially increase max_iterations for better results, but longer runtimes\n", - "warp_matrices, alignment_pairs = imageutils.align_capture(capture,\n", - " ref_index = match_index,\n", - " max_iterations = max_alignment_iterations,\n", - " warp_mode = warp_mode,\n", - " pyramid_levels = pyramid_levels)\n", - "\n", - "print(\"Finished Aligning, warp matrices={}\".format(warp_matrices))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Crop Aligned Images\n", - "After finding image alignments we may need to remove pixels around the edges which aren't present in every image in the capture. To do this we use the affine transforms found above and the image distortions from the image metadata. OpenCV provides a couple of handy helpers for this task in the `cv2.undistortPoints()` and `cv2.transform()` methods. These methods takes a set of pixel coordinates and apply our undistortion matrix and our affine transform, respectively. So, just as we did when registering the images, we first apply the undistortion process the coordinates of the image borders, then we apply the affine transformation to that result. The resulting pixel coordinates tell us where the image borders end up after this pair of transformations, and we can then crop the resultant image to these coordinates." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cropped_dimensions, edges = imageutils.find_crop_bounds(capture, warp_matrices, warp_mode=warp_mode)\n", - "im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, match_index, img_type=img_type)" - ] - }, - { - "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": { - "scrolled": false - }, - "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 = [capture.band_names_lower().index('red'),\n", - " capture.band_names_lower().index('green'),\n", - " capture.band_names_lower().index('blue')]\n", - "cir_band_indices = [capture.band_names_lower().index('nir'),\n", - " capture.band_names_lower().index('red'),\n", - " capture.band_names_lower().index('green')]\n", - "\n", - "# Create a 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[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast\n", - "im_max = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values\n", - "\n", - "# for rgb true color, we use the same min and max scaling across the 3 bands to \n", - "# maintain the \"white balance\" of the calibrated image\n", - "for i in rgb_band_indices:\n", - " im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)\n", - "\n", - "rgb = im_display[:,:,rgb_band_indices]\n", - "\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": {}, - "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": [ - "## Image Export\n", - "\n", - "Composite images can be exported to JPEG or PNG format using the `imageio` package. These images may be useful for visualization or thumbnailing, and creating RGB thumbnails of a set of images can provide a convenient way to browse the imagery in a more visually appealing way that browsing the raw imagery. " - ] - }, - { - "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'))\n", - "imageio.imwrite('cir.'+imtype, (255*cir).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", - "nir_band = capture.band_names_lower().index('nir')\n", - "red_band = capture.band_names_lower().index('red')\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[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])\n", - "\n", - "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", - "if img_type == 'reflectance':\n", - " ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < 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[:,:,nir_band] < 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", - "rededge_band = capture.band_names_lower().index('red edge')\n", - "ndre = (im_aligned[:,:,nir_band] - im_aligned[:,:,rededge_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,rededge_band])\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": {}, - "outputs": [], - "source": [ - "if len(capture.lw_indices()) > 0:\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": [ - "## Red vs NIR Reflectance\n", - "\n", - "Finally, we show a classic agricultural remote sensing output in the tassled cap plot. This plot can be useful for visualizing row crops and plots the Red Reflectance channel on the X-axis against the NIR reflectance channel on the Y-axis. This plot also clearly shows the line of the soil in that space. The tassled cap view isn't very useful for this arid data set; however, we can see the \"badge of trees\" of high NIR reflectance and relatively low red reflectance. This provides an example of one of the uses of aligned images for single capture analysis." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "x_band = red_band\n", - "y_band = nir_band\n", - "x_max = np.max(im_aligned[:,:,x_band])\n", - "y_max = np.max(im_aligned[:,:,y_band])\n", - "\n", - "fig = plt.figure(figsize=(12,12))\n", - "plt.hexbin(im_aligned[:,:,x_band],im_aligned[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))\n", - "ax = fig.gca()\n", - "ax.set_xlim([0,x_max])\n", - "ax.set_ylim([0,y_max])\n", - "plt.xlabel(\"{} Reflectance\".format(capture.band_names()[x_band]))\n", - "plt.ylabel(\"{} Reflectance\".format(capture.band_names()[y_band]))\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Print warp_matrices for usage elsewhere, such as Batch Processing\n", - "Last, we output the warp_matrices that we got for this image stack for usage elsewhere. Currently these can be used in the `Batch Processing.ipynb` notebook to save reflectance-compensated stacks of images to a directory." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "print(warp_matrices)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "---\n", - "Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the [project git repository](https://github.com/micasense/imageprocessing)" - ] - } - ], - "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.7.3" - }, - "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 -} diff --git a/Batch Processing-10Band.ipynb b/Batch Processing-10Band.ipynb deleted file mode 100644 index 4f083b79..00000000 --- a/Batch Processing-10Band.ipynb +++ /dev/null @@ -1,351 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 10-band Batch Processing Example\n", - "\n", - "In this example, we use the `micasense.imageset` class to load a set of directories of images into a list of `micasense.capture` objects, and we iterate over that list saving out each image as an aligned stack of images as separate bands in a single tiff file each. Next, we use the metadata from the original captures to write out a log file of the captures and their locations. Finally, we use `exiftool` from the command line to inject that metadata into the processed images, allowing us to stitch those images using commercial software such as Pix4D or Agisoft.\n", - "\n", - "For an example dataset, download and unzip either sample dataset at https://www.micasense.com/dual-camera-sample-data into your Downloads folder, and ensure the paths below point to the correct location." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load Images into ImageSet" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import FloatProgress, Layout\n", - "from IPython.display import display\n", - "import micasense.imageset as imageset\n", - "import micasense.capture as capture\n", - "import os, glob\n", - "import multiprocessing\n", - "\n", - "panelNames = None\n", - "useDLS = True\n", - "\n", - "imagePath = os.path.expanduser(os.path.join('~','Downloads','DualCam-Farm','farm_only'))\n", - "panelNames = glob.glob(os.path.join(imagePath,'IMG_0002_*.tif'))\n", - "\n", - "outputPath = os.path.join(imagePath,'..','stacks')\n", - "thumbnailPath = os.path.join(outputPath, '..', 'thumbnails')\n", - "\n", - "overwrite = False # Set to False to continue interrupted processing\n", - "generateThumbnails = True\n", - "\n", - "# Allow this code to align both radiance and reflectance images; bu excluding\n", - "# a definition for panelNames above, radiance images will be used\n", - "# For panel images, efforts will be made to automatically extract the panel information\n", - "# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance\n", - "# will need to be set in the panel_reflectance_by_band variable.\n", - "# Note: radiance images will not be used to properly create NDVI/NDRE images below.\n", - "if panelNames is not None:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "else:\n", - " panelCap = None\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.65]*len(panelCap.images) #inexact, but quick\n", - " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band) \n", - " img_type = \"reflectance\"\n", - "else:\n", - " if useDLS:\n", - " img_type='reflectance'\n", - " else:\n", - " img_type = \"radiance\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## This progress widget is used for display of the long-running process\n", - "f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Loading\")\n", - "display(f)\n", - "def update_f(val):\n", - " if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x\n", - " f.value=val\n", - "\n", - "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "update_f(1.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import math\n", - "import numpy as np\n", - "from mapboxgl.viz import *\n", - "from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between\n", - "from mapboxgl.utils import create_color_stops\n", - "import pandas as pd\n", - "\n", - "data, columns = imgset.as_nested_lists()\n", - "df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)\n", - "\n", - "#Insert your mapbox token here\n", - "token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'\n", - "color_property = 'dls-yaw'\n", - "color_property = 'altitude'\n", - "num_color_classes = 8\n", - "\n", - "min_val = df[color_property].min()\n", - "max_val = df[color_property].max()\n", - "\n", - "import jenkspy\n", - "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", - "\n", - "color_stops = create_color_stops(breaks,colors='YlOrRd')\n", - "geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')\n", - "\n", - "viz = CircleViz(geojson_data, access_token=token, color_property=color_property,\n", - " color_stops=color_stops,\n", - " center=[df['longitude'].median(),df['latitude'].median()], \n", - " zoom=16, height='600px',\n", - " style='mapbox://styles/mapbox/satellite-streets-v9')\n", - "viz.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Define which warp method to use\n", - "For newer data sets with RigRelatives tags (images captured with RedEdge version 3.4.0 or greater with a valid calibration load, see https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields), we can use the RigRelatives for a simple alignment.\n", - "\n", - "For sets without those tags, or sets that require a RigRelatives optimization, we can go through the Alignment.ipynb notebook and get a set of `warp_matrices` that we can use here to align." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from numpy import array\n", - "from numpy import float32\n", - "\n", - "# Use the warp_matrices derived from the Alignment Tutorial for this RedEdge set without RigRelatives\n", - "warp_matrices = [array([[ 1.0020243e+00, -3.7388311e-04, 2.4971788e+01],\n", - " [ 6.7297497e-04, 1.0005866e+00, 1.7188536e+01],\n", - " [ 2.4259109e-06, -9.2373267e-07, 1.0000000e+00]], dtype=float32), array([[ 9.9140632e-01, -4.6332614e-05, 4.8500401e+01],\n", - " [ 3.2340995e-05, 9.9200422e-01, -1.0915921e+01],\n", - " [-7.3704086e-07, 5.0890253e-07, 1.0000000e+00]], dtype=float32), array([[ 1.0018263e+00, -2.1731904e-04, 5.5316315e+00],\n", - " [ 7.2411756e-04, 1.0021795e+00, 5.8745198e+00],\n", - " [-1.9047379e-08, 9.7758209e-07, 1.0000000e+00]], dtype=float32), array([[ 9.9152303e-01, -5.4825414e-03, 4.1536880e+01],\n", - " [ 3.8441001e-03, 9.9495757e-01, 1.7250452e+01],\n", - " [-3.2921032e-06, -2.4233820e-08, 1.0000000e+00]], dtype=float32), array([[ 1.0006192e+00, -3.0658240e-04, -2.5816131e-01],\n", - " [ 7.8755329e-05, 9.9954307e-01, 2.9809377e-01],\n", - " [ 9.1640561e-07, -1.0784843e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9773926e-01, -6.3800282e-04, 5.2199936e+01],\n", - " [-3.4246168e-03, 9.9601907e-01, 2.0550659e+01],\n", - " [-4.6251063e-07, -4.8716843e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9622118e-01, 3.1637053e-03, 3.7498917e+01],\n", - " [-6.7951437e-03, 9.9743211e-01, 8.9517927e+00],\n", - " [-3.6472218e-06, -2.4649705e-06, 1.0000000e+00]], dtype=float32), array([[ 9.8943901e-01, 3.7658634e-04, 9.4948044e+00],\n", - " [-4.0384033e-03, 9.8851675e-01, 1.5366467e+01],\n", - " [-2.4371677e-06, -3.8438825e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9749213e-01, 1.6272087e-03, 4.3243721e-01],\n", - " [-7.3282972e-05, 9.9533182e-01, 3.5523354e+01],\n", - " [ 3.8597086e-06, -4.0187538e-07, 1.0000000e+00]], dtype=float32), array([[ 9.9992698e-01, 6.6664284e-03, -9.0784521e+00],\n", - " [-9.0053231e-03, 9.9836856e-01, 1.5190173e+01],\n", - " [-1.6761204e-07, -3.6131762e-06, 1.0000000e+00]], dtype=float32)]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Align images and save each capture to a layered tiff file" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import exiftool\n", - "import datetime\n", - "\n", - "use_multi_process = True # set to False for single-process saving\n", - "overwrite_existing = False # skip existing files, set to True to overwrite\n", - "\n", - "## This progress widget is used for display of the long-running process\n", - "f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Saving\")\n", - "display(f2)\n", - "def update_f2(val):\n", - " f2.value=val\n", - "\n", - "if not os.path.exists(outputPath):\n", - " os.makedirs(outputPath)\n", - "if generateThumbnails and not os.path.exists(thumbnailPath):\n", - " os.makedirs(thumbnailPath)\n", - "\n", - "# Save out geojson data so we can open the image capture locations in our GIS\n", - "with open(os.path.join(outputPath,'imageSet.json'),'w') as f:\n", - " f.write(str(geojson_data))\n", - "\n", - "# If we didn't provide a panel above, irradiance set to None will cause DLS data to be used\n", - "try:\n", - " irradiance = panel_irradiance+[0]\n", - "except NameError:\n", - " irradiance = None\n", - "\n", - "start_time = datetime.datetime.now()\n", - "\n", - "# Save all captures in the imageset as aligned stacks\n", - "imgset.save_stacks(warp_matrices,\n", - " outputPath,\n", - " thumbnailPath,\n", - " irradiance = irradiance,\n", - " multiprocess=use_multi_process, \n", - " overwrite=overwrite_existing, \n", - " progress_callback=update_f2)\n", - "\n", - "end_time = datetime.datetime.now()\n", - "update_f2(1.0)\n", - "\n", - "print(\"Saving time: {}\".format(end_time-start_time))\n", - "print(\"Alignment+Saving rate: {:.2f} captures per second\".format(float(len(imgset.captures))/float((end_time-start_time).total_seconds())))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Extract Metadata from Captures list and save to log.csv" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def decdeg2dms(dd):\n", - " is_positive = dd >= 0\n", - " dd = abs(dd)\n", - " minutes,seconds = divmod(dd*3600,60)\n", - " degrees,minutes = divmod(minutes,60)\n", - " degrees = degrees if is_positive else -degrees\n", - " return (degrees,minutes,seconds)\n", - "\n", - "header = \"SourceFile,\\\n", - "GPSDateStamp,GPSTimeStamp,\\\n", - "GPSLatitude,GpsLatitudeRef,\\\n", - "GPSLongitude,GPSLongitudeRef,\\\n", - "GPSAltitude,GPSAltitudeRef,\\\n", - "FocalLength,\\\n", - "XResolution,YResolution,ResolutionUnits\\n\"\n", - "\n", - "lines = [header]\n", - "for capture in imgset.captures:\n", - " #get lat,lon,alt,time\n", - " outputFilename = capture.uuid+'.tif'\n", - " fullOutputPath = os.path.join(outputPath, outputFilename)\n", - " lat,lon,alt = capture.location()\n", - " #write to csv in format:\n", - " # IMG_0199_1.tif,\"33 deg 32' 9.73\"\" N\",\"111 deg 51' 1.41\"\" W\",526 m Above Sea Level\n", - " latdeg, latmin, latsec = decdeg2dms(lat)\n", - " londeg, lonmin, lonsec = decdeg2dms(lon)\n", - " latdir = 'North'\n", - " if latdeg < 0:\n", - " latdeg = -latdeg\n", - " latdir = 'South'\n", - " londir = 'East'\n", - " if londeg < 0:\n", - " londeg = -londeg\n", - " londir = 'West'\n", - " resolution = capture.images[0].focal_plane_resolution_px_per_mm\n", - "\n", - " linestr = '\"{}\",'.format(fullOutputPath)\n", - " linestr += capture.utc_time().strftime(\"%Y:%m:%d,%H:%M:%S,\")\n", - " linestr += '\"{:d} deg {:d}\\' {:.2f}\"\" {}\",{},'.format(int(latdeg),int(latmin),latsec,latdir[0],latdir)\n", - " linestr += '\"{:d} deg {:d}\\' {:.2f}\"\" {}\",{},{:.1f} m Above Sea Level,Above Sea Level,'.format(int(londeg),int(lonmin),lonsec,londir[0],londir,alt)\n", - " linestr += '{}'.format(capture.images[0].focal_length)\n", - " linestr += '{},{},mm'.format(resolution,resolution)\n", - " linestr += '\\n' # when writing in text mode, the write command will convert to os.linesep\n", - " lines.append(linestr)\n", - "\n", - "fullCsvPath = os.path.join(outputPath,'log.csv')\n", - "with open(fullCsvPath, 'w') as csvfile: #create CSV\n", - " csvfile.writelines(lines)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Use Exiftool from the command line to write metadata to images" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import subprocess\n", - "\n", - "if os.environ.get('exiftoolpath') is not None:\n", - " exiftool_cmd = os.path.normpath(os.environ.get('exiftoolpath'))\n", - "else:\n", - " exiftool_cmd = 'exiftool'\n", - " \n", - "cmd = '{} -csv=\"{}\" -overwrite_original {}'.format(exiftool_cmd, fullCsvPath, outputPath)\n", - "print(cmd)\n", - "if(subprocess.check_call(cmd) == 0):\n", - " print(\"Successfully updated stack metadata\")" - ] - } - ], - "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.7.6" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/Batch Processing.ipynb b/Batch Processing.ipynb deleted file mode 100644 index 2c5557e7..00000000 --- a/Batch Processing.ipynb +++ /dev/null @@ -1,341 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Batch Processing Example\n", - "\n", - "In this example, we use the `micasense.imageset` class to load a set of directories of images into a list of `micasense.capture` objects, and we iterate over that list saving out each image as an aligned stack of images as separate bands in a single tiff file each. Next, we use the metadata from the original captures to write out a log file of the captures and their locations. Finally, we use `exiftool` from the command line to inject that metadata into the processed images, allowing us to stitch those images using commercial software such as Pix4D or Agisoft.\n", - "\n", - "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge version of at least 3.4.0 or any version of Altum. If your images don't meet that spec, you can also follow this support ticket to add the RigRelatives tag to them: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load Images into ImageSet" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import FloatProgress, Layout\n", - "from IPython.display import display\n", - "import micasense.imageset as imageset\n", - "import micasense.capture as capture\n", - "import os, glob\n", - "import multiprocessing\n", - "\n", - "panelNames = None\n", - "useDLS = True\n", - "\n", - "imagePath = os.path.expanduser(os.path.join('~','Downloads','RedEdgeImageSet','0000SET'))\n", - "panelNames = glob.glob(os.path.join(imagePath,'000','IMG_0000_*.tif'))\n", - "panelCap = capture.Capture.from_filelist(panelNames)\n", - "\n", - "outputPath = os.path.join(imagePath,'..','stacks')\n", - "thumbnailPath = os.path.join(outputPath, '..', 'thumbnails')\n", - "\n", - "overwrite = False # can be set to set to False to continue interrupted processing\n", - "generateThumbnails = True\n", - "\n", - "# Allow this code to align both radiance and reflectance images; bu excluding\n", - "# a definition for panelNames above, radiance images will be used\n", - "# For panel images, efforts will be made to automatically extract the panel information\n", - "# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance\n", - "# will need to be set in the panel_reflectance_by_band variable.\n", - "# Note: radiance images will not be used to properly create NDVI/NDRE images below.\n", - "if panelNames is not None:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "else:\n", - " panelCap = None\n", - "\n", - "if panelCap is not None:\n", - " if panelCap.panel_albedo() is not None and not any(v is None for v in panelCap.panel_albedo()):\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] #RedEdge band_index order\n", - " \n", - " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band) \n", - " img_type = \"reflectance\"\n", - "else:\n", - " if useDLS:\n", - " img_type='reflectance'\n", - " else:\n", - " img_type = \"radiance\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## This progress widget is used for display of the long-running process\n", - "f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Loading\")\n", - "display(f)\n", - "def update_f(val):\n", - " if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x\n", - " f.value=val\n", - "\n", - "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "update_f(1.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import math\n", - "import numpy as np\n", - "from mapboxgl.viz import *\n", - "from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between\n", - "from mapboxgl.utils import create_color_stops\n", - "import pandas as pd\n", - "\n", - "data, columns = imgset.as_nested_lists()\n", - "df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)\n", - "\n", - "#Insert your mapbox token here\n", - "token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'\n", - "color_property = 'dls-yaw'\n", - "num_color_classes = 8\n", - "\n", - "min_val = df[color_property].min()\n", - "max_val = df[color_property].max()\n", - "\n", - "import jenkspy\n", - "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", - "\n", - "color_stops = create_color_stops(breaks,colors='YlOrRd')\n", - "geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')\n", - "\n", - "viz = CircleViz(geojson_data, access_token=token, color_property=color_property,\n", - " color_stops=color_stops,\n", - " center=[df['longitude'].median(),df['latitude'].median()], \n", - " zoom=16, height='600px',\n", - " style='mapbox://styles/mapbox/satellite-streets-v9')\n", - "viz.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Define which warp method to use\n", - "For newer data sets with RigRelatives tags (images captured with RedEdge version 3.4.0 or greater with a valid calibration load, see https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields), we can use the RigRelatives for a simple alignment.\n", - "\n", - "For sets without those tags, or sets that require a RigRelatives optimization, we can go through the Alignment.ipynb notebook and get a set of `warp_matrices` that we can use here to align." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from numpy import array\n", - "from numpy import float32\n", - "\n", - "# Set warp_matrices to none to align using RigRelatives\n", - "# Or\n", - "# Use the warp_matrices derived from the Alignment Tutorial for this RedEdge set without RigRelatives\n", - "warp_matrices = [array([[ 1.0022864e+00, -2.5218755e-03, -7.8898020e+00],\n", - " [ 2.3614739e-03, 1.0036649e+00, -1.3134377e+01],\n", - " [-1.7785899e-06, 1.1343118e-06, 1.0000000e+00]], dtype=float32), array([[1., 0., 0.],\n", - " [0., 1., 0.],\n", - " [0., 0., 1.]], dtype=float32), array([[ 9.9724638e-01, -1.5535230e-03, 1.2301294e+00],\n", - " [ 8.6745428e-04, 9.9738181e-01, -1.6499169e+00],\n", - " [-8.2816513e-07, -3.4488804e-07, 1.0000000e+00]], dtype=float32), array([[ 1.0007139e+00, -8.4427800e-03, 1.6312805e+01],\n", - " [ 6.2834378e-03, 9.9977130e-01, -1.6011697e+00],\n", - " [-1.9520389e-06, -6.3762940e-07, 1.0000000e+00]], dtype=float32), array([[ 9.9284178e-01, 9.2155562e-04, 1.6069822e+01],\n", - " [-3.2895457e-03, 9.9262553e-01, -5.0333548e-01],\n", - " [-1.5845577e-06, -1.7680986e-06, 1.0000000e+00]], dtype=float32)]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Align images and save each capture to a layered tiff file" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import exiftool\n", - "import datetime\n", - "## This progress widget is used for display of the long-running process\n", - "f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Saving\")\n", - "display(f2)\n", - "def update_f2(val):\n", - " f2.value=val\n", - "\n", - "if not os.path.exists(outputPath):\n", - " os.makedirs(outputPath)\n", - "if generateThumbnails and not os.path.exists(thumbnailPath):\n", - " os.makedirs(thumbnailPath)\n", - "\n", - "# Save out geojson data so we can open the image capture locations in our GIS\n", - "with open(os.path.join(outputPath,'imageSet.json'),'w') as f:\n", - " f.write(str(geojson_data))\n", - " \n", - "try:\n", - " irradiance = panel_irradiance+[0]\n", - "except NameError:\n", - " irradiance = None\n", - "\n", - "start = datetime.datetime.now()\n", - "for i,capture in enumerate(imgset.captures):\n", - " outputFilename = capture.uuid+'.tif'\n", - " thumbnailFilename = capture.uuid+'.jpg'\n", - " fullOutputPath = os.path.join(outputPath, outputFilename)\n", - " fullThumbnailPath= os.path.join(thumbnailPath, thumbnailFilename)\n", - " if (not os.path.exists(fullOutputPath)) or overwrite:\n", - " if(len(capture.images) == len(imgset.captures[0].images)):\n", - " capture.create_aligned_capture(irradiance_list=irradiance, warp_matrices=warp_matrices)\n", - " capture.save_capture_as_stack(fullOutputPath)\n", - " if generateThumbnails:\n", - " capture.save_capture_as_rgb(fullThumbnailPath)\n", - " capture.clear_image_data()\n", - " update_f2(float(i)/float(len(imgset.captures)))\n", - "update_f2(1.0)\n", - "end = datetime.datetime.now()\n", - "\n", - "print(\"Saving time: {}\".format(end-start))\n", - "print(\"Alignment+Saving rate: {:.2f} images per second\".format(float(len(imgset.captures))/float((end-start).total_seconds())))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Extract Metadata from Captures list and save to log.csv" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def decdeg2dms(dd):\n", - " is_positive = dd >= 0\n", - " dd = abs(dd)\n", - " minutes,seconds = divmod(dd*3600,60)\n", - " degrees,minutes = divmod(minutes,60)\n", - " degrees = degrees if is_positive else -degrees\n", - " return (degrees,minutes,seconds)\n", - "\n", - "header = \"SourceFile,\\\n", - "GPSDateStamp,GPSTimeStamp,\\\n", - "GPSLatitude,GpsLatitudeRef,\\\n", - "GPSLongitude,GPSLongitudeRef,\\\n", - "GPSAltitude,GPSAltitudeRef,\\\n", - "FocalLength,\\\n", - "XResolution,YResolution,ResolutionUnits\\n\"\n", - "\n", - "lines = [header]\n", - "for capture in imgset.captures:\n", - " #get lat,lon,alt,time\n", - " outputFilename = capture.uuid+'.tif'\n", - " fullOutputPath = os.path.join(outputPath, outputFilename)\n", - " lat,lon,alt = capture.location()\n", - " #write to csv in format:\n", - " # IMG_0199_1.tif,\"33 deg 32' 9.73\"\" N\",\"111 deg 51' 1.41\"\" W\",526 m Above Sea Level\n", - " latdeg, latmin, latsec = decdeg2dms(lat)\n", - " londeg, lonmin, lonsec = decdeg2dms(lon)\n", - " latdir = 'North'\n", - " if latdeg < 0:\n", - " latdeg = -latdeg\n", - " latdir = 'South'\n", - " londir = 'East'\n", - " if londeg < 0:\n", - " londeg = -londeg\n", - " londir = 'West'\n", - " resolution = capture.images[0].focal_plane_resolution_px_per_mm\n", - "\n", - " linestr = '\"{}\",'.format(fullOutputPath)\n", - " linestr += capture.utc_time().strftime(\"%Y:%m:%d,%H:%M:%S,\")\n", - " linestr += '\"{:d} deg {:d}\\' {:.2f}\"\" {}\",{},'.format(int(latdeg),int(latmin),latsec,latdir[0],latdir)\n", - " linestr += '\"{:d} deg {:d}\\' {:.2f}\"\" {}\",{},{:.1f} m Above Sea Level,Above Sea Level,'.format(int(londeg),int(lonmin),lonsec,londir[0],londir,alt)\n", - " linestr += '{}'.format(capture.images[0].focal_length)\n", - " linestr += '{},{},mm'.format(resolution,resolution)\n", - " linestr += '\\n' # when writing in text mode, the write command will convert to os.linesep\n", - " lines.append(linestr)\n", - "\n", - "fullCsvPath = os.path.join(outputPath,'log.csv')\n", - "with open(fullCsvPath, 'w') as csvfile: #create CSV\n", - " csvfile.writelines(lines)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Use Exiftool from the command line to write metadata to images" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import subprocess\n", - "\n", - "if os.environ.get('exiftoolpath') is not None:\n", - " exiftool_cmd = os.path.normpath(os.environ.get('exiftoolpath'))\n", - "else:\n", - " exiftool_cmd = 'exiftool'\n", - " \n", - "cmd = '{} -csv=\"{}\" -overwrite_original {}'.format(exiftool_cmd, fullCsvPath, outputPath)\n", - "print(cmd)\n", - "subprocess.check_call(cmd)" - ] - } - ], - "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.7.3" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/Captures.ipynb b/Captures.ipynb deleted file mode 100644 index 9b02d519..00000000 --- a/Captures.ipynb +++ /dev/null @@ -1,117 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Captures\n", - "A `micasense.capture.Capture` object holds a set of 5 (or in the case of Altum 6) images which are captured at the same moment together in a MicaSense camera. Files which meet this criteria will tend to have the same filename except for the suffix, but that is not required to load a captures. Captures can be loaded by starting with one image by calling Capture.from_file('file_name.tif') and adding others using the append_file, or by providing a list of filenames or images. See capture.py for more creation methods." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os, glob\n", - "import micasense.capture as capture\n", - "%matplotlib inline\n", - "\n", - "imagePath = os.path.join('.','data','0000SET','000')\n", - "imageNames = glob.glob(os.path.join(imagePath,'IMG_0000_*.tif'))\n", - "\n", - "capture = capture.Capture.from_filelist(imageNames)\n", - "capture.plot_raw()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# More Capture visualization functions" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "capture.plot_vignette();\n", - "capture.plot_undistorted_radiance();\n", - "capture.plot_panels()" - ] - }, - { - "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", - "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.1" - }, - "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 -} diff --git a/ImageSets.ipynb b/ImageSets.ipynb deleted file mode 100644 index f0b9d4af..00000000 --- a/ImageSets.ipynb +++ /dev/null @@ -1,303 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# ImageSet Extended Examples\n", - "\n", - "A large group of images captured using a RedEdge over a central California orchard are available for [download here](https://s3-us-west-2.amazonaws.com/sample.micasense.com/imageprocessing/RedEdgeImageSet.zip).\n", - "\n", - "With this set extracted to a working folder, this notebook will load all of the images in the set and provide more usage examples of ImageSet data.\n", - "\n", - "## Load a directory of images into an ImageSet" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "from ipywidgets import FloatProgress, Layout\n", - "from IPython.display import display\n", - "import micasense.imageset as imageset\n", - "import os\n", - "\n", - "## This progress widget is used for display of the long-running process\n", - "f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Loading\")\n", - "display(f)\n", - "def update_f(val):\n", - " if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x\n", - " f.value=val\n", - "\n", - "images_dir = os.path.expanduser(os.path.join('~','Downloads','RedEdgeImageSet','0000SET'))\n", - "%time imgset = imageset.ImageSet.from_directory(images_dir, progress_callback=update_f)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Import Capture data into Pandas DataFrame" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import pandas as pd\n", - "data, columns = imgset.as_nested_lists()\n", - "print(\"Columns: {}\".format(columns))\n", - "df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Mapping capture locations using MapboxGL\n", - "\n", - "Below we use the `mapboxgl` extension to plot the measured DLS yaw (or heading) angle from each image's meatadata over the whole imageset. We draw circles at each image location, and then color the circle based on the yaw value. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import math\n", - "import numpy as np\n", - "from mapboxgl.viz import *\n", - "from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between\n", - "from mapboxgl.utils import create_color_stops\n", - "\n", - "#Insert your mapbox token here\n", - "token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'\n", - "color_stops = create_color_stops(np.linspace(-math.pi,math.pi,num=8),colors='YlOrRd')\n", - "data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')\n", - "viz = CircleViz(data, access_token=token, color_property='dls-yaw',\n", - " color_stops=color_stops,\n", - " center=[df['longitude'].median(),df['latitude'].median()], \n", - " zoom=16, height='600px',\n", - " style='mapbox://styles/mapbox/satellite-streets-v9')\n", - "viz.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plotting capture metadata" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "%matplotlib inline\n", - "ax=df.plot(y=columns[3:], subplots=True, figsize=(16, 10), style=['g','c','y','k','b','g','r','m','k']);\n", - "for a in ax:\n", - " a.legend(loc='right', bbox_to_anchor=(1.1, 0.5), ncol=1, fancybox=True, shadow=True);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Separating ground and flight images\n", - "\n", - "We can plot a histogram of the image altitudes for all images in the flight and use the basic statistics to guess at an altitude below which the panel images were likely taken. This can give us an altitude threshold we can apply to separate images taken in flight and on the ground." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot the histogram of the altitude data\n", - "df.altitude.hist();\n", - "# find the altitude above which the flight images occur\n", - "cutoff_altitude = df.altitude.mean()-3.0*df.altitude.std()\n", - "plt.axvline(x=cutoff_altitude,c='r');\n", - "plt.xlabel('Capture altitude (m)');\n", - "plt.ylabel('Number of occurances');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Automatically extracting radiances of panel images\n", - "Using the Panel class, we can automatically find and compute the radiance of panel regions within panel images. Each Capture object exposes a `panel_radiance()` method which can be used on panel images to get a list of the radiance of each band (in the native RedEdge band order)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "flight = df.altitude>cutoff_altitude\n", - "ground = ~flight\n", - "ground_idx = np.arange(len(ground))[ground]\n", - "flight_idx = np.arange(len(ground))[flight]\n", - "ground_captures = np.array(imgset.captures)[ground_idx]\n", - "flight_captures = np.array(imgset.captures)[flight_idx]\n", - "\n", - "panel_radiances = []\n", - "dls_irradiances = []\n", - "panel_timestamps = []\n", - "for cap in ground_captures:\n", - " if cap.panels_in_all_expected_images():\n", - " panel_timestamps.append(cap.utc_time())\n", - " panel_radiances.append(cap.panel_radiance())\n", - " dls_irradiances.append(cap.dls_irradiance())\n", - "\n", - "dls_irradiances = np.asarray(dls_irradiances)\n", - "panel_radiances = np.asarray(panel_radiances)\n", - "\n", - "###\n", - "panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] #RedEdge band_index order\n", - "panel_irradiance = ground_captures[0].panel_irradiance(panel_reflectance_by_band)\n", - "plt.figure()\n", - "plt.scatter(ground_captures[0].center_wavelengths(), panel_irradiance);\n", - "plt.xscale('log')\n", - "plt.ylabel(\"Irradiance (w/m^2/nm)\");\n", - "plt.xlabel(\"Wavelength (nm)\");\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plotting DLS Irradiances\n", - "Once we have extracted the panel radiances from the ground images, we can convert them to irradiance as we did in tutorial #1 and plot the irradiances over time and extract radiance to reflectance conversion factors." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import math\n", - "\n", - "df[df.altitude>cutoff_altitude].plot(y=columns[8:13], figsize=(14,8))\n", - "plt.ylabel(\"Irradiance (w/m^2/nm)\");\n" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "collapsed": true - }, - "source": [ - "---\n", - "Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the [project git repository](https://github.com/micasense/imageprocessing)" - ] - } - ], - "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.7.1" - }, - "toc": { - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "toc_cell": false, - "toc_position": { - "height": "1238px", - "left": "0px", - "right": "1114.56px", - "top": "111px", - "width": "383px" - }, - "toc_section_display": "block", - "toc_window_display": true - }, - "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()) " - } - }, - "position": { - "height": "144px", - "left": "1129px", - "right": "20px", - "top": "119px", - "width": "350px" - }, - "types_to_exclude": [ - "module", - "function", - "builtin_function_or_method", - "instance", - "_Feature" - ], - "window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/Images.ipynb b/Images.ipynb deleted file mode 100644 index efb29c99..00000000 --- a/Images.ipynb +++ /dev/null @@ -1,145 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "collapsed": true - }, - "source": [ - "# Images\n", - "This notebook shows usage for the Image type. This type is useful for loading single image files, reading the embedded metadata, and performing radiometric calibrations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from micasense.image import Image\n", - "import os, glob\n", - "%matplotlib inline\n", - "\n", - "imagePath = os.path.join('.','data','0000SET','000')\n", - "imageName = glob.glob(os.path.join(imagePath,'IMG_0000_1.tif'))[0]\n", - "\n", - "img = Image(imageName)\n", - "img.plot_all();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Image Histogram" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import cv2\n", - "\n", - "nbins = 1024\n", - "vmin = 0\n", - "vmax = 2**16\n", - "bins = range(vmin,vmax, int(vmax/nbins))\n", - "hist = cv2.calcHist([img.raw().ravel()],[0],None,[nbins],[vmin,vmax])\n", - "plt.plot(bins,hist);\n", - "plt.xlim(img.raw().min(),img.raw().max())\n", - "plt.ylim(0,hist.max())\n", - "plt.xlabel('Pixel Value')\n", - "plt.ylabel('Frequency')\n", - "plt.show()\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Captures\n", - "\n", - "Once you've mastered images, move on to groups of images, or [Captures](Captures.html)" - ] - }, - { - "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", - "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.1" - }, - "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()) " - } - }, - "position": { - "height": "285px", - "left": "453px", - "right": "20px", - "top": "105px", - "width": "451px" - }, - "types_to_exclude": [ - "module", - "function", - "builtin_function_or_method", - "instance", - "_Feature" - ], - "window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/MicaSense Image Processing Setup.ipynb b/MicaSense Image Processing Setup.ipynb deleted file mode 100755 index bdf8d3ea..00000000 --- a/MicaSense Image Processing Setup.ipynb +++ /dev/null @@ -1,288 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "tags": [] - }, - "source": [ - "# MicaSense Image Processing Setup\n", - "\n", - "## Overview\n", - "\n", - "This series of tutorials will be a walk through on how to process RedEdge data from raw images through conversion to reflectance. In this first tutorial, we will cover the tools required to do this, get them installed, and verify that the installation works.\n", - "\n", - "## System Requirements\n", - "\n", - "Our tutorials are written using Python3. Python has great library support for image processing through libraries such as OpenCV, SciKit Image, and others. In this tutorial, we'll use Python, OpenCV, numpy, and matplotlib, as well as the standalone exiftool and it's Python wrapper to open and manipulate RedEdge images to transform raw digital number values into quantitative reflectance. Python 2.7 can work for this tutorial, but we only test with Python3 and later tutorials use some extra libraries that are best supported in Python3 (specifically pysolar), so we recommend that if you're starting with Python from scratch to install Python3.\n", - "\n", - "This tutorial has been tested on Windows, MacOS, and Linux. It is likely to work on other platforms, especially unix-based platforms like macOS, but you will have to do the legwork to get the required software installed and working.\n", - "\n", - "### Software/Libraries Overview\n", - "\n", - "The following softare and libraries are required for this tutorial:\n", - "\n", - "* [python3](https://www.python.org/downloads/release/latest)\n", - "* [numpy](https://www.numpy.org/)\n", - "* [openCV](https://opencv.org/releases.html)\n", - "* [matplotlib](https://matplotlib.org/users/installing.html)\n", - "* [exiftool](https://exiftool.org/) + [pyexiftool](https://github.com/sylikc/pyexiftool)\n", - "* [scikit-image](https://scikit-image.org/)\n", - "* [zbar](http://zbar.sourceforge.net/) + [pyzbar](https://github.com/NaturalHistoryMuseum/pyzbar)\n", - "* [pysolar](http://pysolar.org/)\n", - "* [pandas](https://pandas.pydata.org/)\n", - "* [mapboxgl](https://github.com/mapbox/mapboxgl-jupyter)\n", - "\n", - "Below, we go through the options to download and install a full working Python environment with these tools (and their dependencies). We're using the [Anaconda](https://www.anaconda.com/download/) or [miniconda](https://conda.io/miniconda.html) environments where possible to ease installation, but if you're already a Python package management guru, you can use `git` to checkout this code repository and look at the `micasense_conda_env.yml` file for the dependencies you'll need in your virtual environment.\n", - "\n", - "### Linux (Debian/Ubuntu)\n", - "\n", - "For linux (and Mac, to some extent) you can either install the libraries directly using `pip` or install `miniconda` or `anaconda` to create completely separate environments. We have had success installing `miniconda` locally -- it's a smaller install than `anaconda` and can be installed without using `sudo` and doesn't impact the system-installed Python or Python libraries. You will likely still need to use `sudo` to install \n", - "\n", - "The following is what we had to do on a fresh Ubuntu 18.04 image to install the library. First we installed some system tools and libraries:\n", - "\n", - " sudo apt install git\n", - " sudo apt install libzbar0\n", - " sudo apt install make\n", - " \n", - "Next we installed [exiftool](https://exiftool.org/):\n", - "\n", - " wget https://cpan.metacpan.org/authors/id/E/EX/EXIFTOOL/Image-ExifTool-12.15.tar.gz\n", - " tar -xvzf Image-ExifTool-12.15.tar.gz \n", - " cd Image-ExifTool-12.15/\n", - " perl Makefile.PL \n", - " make test\n", - " sudo make install\n", - "\n", - "Then we installed miniconda. Navigate to the [miniconda download page](https://conda.io/miniconda.html) and download the installer for your system and follow the [installation instructions](https://conda.io/docs/user-guide/install/index.html)\n", - "\n", - "Once these tools are installed, you can check out this repository and create the `micasense conda` environment:\n", - "\n", - " git clone https://github.com/micasense/imageprocessing.git\n", - " cd imageprocessing\n", - " conda env create -f micasense_conda_env.yml\n", - "\n", - "Finally, one way to verify our install by running the built in tests:\n", - "\n", - " cd imageprocessing\n", - " conda activate micasense\n", - " pytest .\n", - "\n", - "Or, to start working with the notebooks (including running the test code below):\n", - "\n", - " cd imageprocessing\n", - " conda activate micasense\n", - " jupyter notebook .\n", - "\n", - "\n", - "### Windows setup\n", - "\n", - "When installing on Windows we rely on the [Anaconda](https://www.anaconda.com/download/) Python environment to do most of the heavy lifting for us.\n", - "\n", - "* Install [Anaconda](https://www.anaconda.com/download/) for your system by downloading the **Python 3.6** version\n", - "\n", - " * When installing Anaconda, choose **\"install for only me\"** instead of \"install for all users,\" as this simplifies installation of other packages\n", - "\n", - "* Download the [exiftool windows package](https://exiftool.org/) and unzip it to a permanent location such as `c:\\exiftool\\`. Now we need to tell the Python code where to find exiftool (so we don't have to set it up in every script we write), and we do that by adding the path to exiftool as an environment variable.\n", - " * Create an environment variable called `exiftoolpath` with a value of the full path to exiftool. For example, `c:\\exiftool\\exiftool.exe`\n", - " * To do this on Windows 10, press Start button or the Windows key, then type `Path` and click `Edit Environment Variables for Your Account`\n", - " * Click `New`\n", - " * In Variable Name type `exiftoolpath`\n", - " * In Variable Value type `c:\\exiftool\\exiftool.exe`\n", - "\n", - "* Open an Anaconda console from the start menu as an administrator by clicking `Start->Anaconda`, right-click `Anaconda Console`, choose `Run as Administrator`. Execute the following commands in the anaconda console:\n", - "\n", - " * `cd` to the directory you git cloned this repository to\n", - " * `conda env create -f micasense_conda_env.yml`\n", - " * This will configure an anaconda environment with all of the required tools and libraries This will take a while depending on your computer and internet speeds (5-10 minutes is not uncommon)\n", - " * When it's done, run `activate micasense` to activate the environment configured\n", - " * Each time you run start a new anaconda prompt, you'll need to run `activate micasense`\n", - " \n", - "### macOS setup\n", - "\n", - "First, install [Homebrew](https://brew.sh/) - The Missing Package Manager for macOS (or Linux). If not already installed, Homebrew will also install macOS `xcode` developer tools for you. Open your Terminal and:\n", - "\n", - "`/bin/bash -c \"$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)\"`\n", - "\n", - "Next, using Homebrew, install the dependencies like:\n", - "\n", - "```\n", - "brew install cask\n", - "brew install exiftool\n", - "brew install zbar\n", - "brew install miniconda # or brew install anaconda \n", - "brew install git\n", - "brew install git-lfs\n", - "```\n", - "\n", - "If instead you're more comfortable with graphical installers, the [Anaconda](https://www.anaconda.com/download/) version for **Python 3.7** may be right for you.\n", - "\n", - "Once these tools are installed, you can check out this repository and create the `micasense conda` environment by opening a Terminal and running the following commands:\n", - " \n", - " git clone https://github.com/micasense/imageprocessing.git\n", - " cd imageprocessing\n", - " conda env create -f micasense_conda_env.yml\n", - " \n", - "This will take a while (5-10 minutes isn't uncommon). Once it's done, one way to verify our install by running the built-in tests:\n", - "\n", - " conda activate micasense\n", - " pytest .\n", - "\n", - "Or, to start working with the notebooks (including running the test code below):\n", - "\n", - " cd imageprocessing\n", - " conda activate micasense\n", - " jupyter notebook .\n", - "\n", - "\n", - "## Running the notebooks\n", - "\n", - "* If running on Windows, run the `Anaconda Prompt` from the Start menu and type `activate micasense`\n", - "* `cd` to the imageprocessing checkout directory\n", - "* Run `jupyter notebook .`\n", - "\n", - "## Testing Installation\n", - "\n", - "The following Python snippet can be run from a jupyter notebook, inside iPython, or by saving to a script and running from the command line. If you're on windows, make sure you have set the location of exiftool in the `exiftoolpath` environment variable. If this script succeeds, your system is ready to go! If not, check the installation documentation for the module import that is having issues.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "import cv2 #openCV\n", - "import exiftool\n", - "import os, glob\n", - "import numpy as np\n", - "import pyzbar.pyzbar as pyzbar\n", - "import matplotlib.pyplot as plt\n", - "import mapboxgl\n", - "\n", - "print()\n", - "print(\"Successfully imported all required libraries.\")\n", - "print()\n", - "\n", - "if os.name == 'nt':\n", - " if os.environ.get('exiftoolpath') is None:\n", - " print(\"Set the `exiftoolpath` environment variable as described above\")\n", - " else:\n", - " if not os.path.isfile(os.environ.get('exiftoolpath')):\n", - " print(\"The provided exiftoolpath isn't a file, check the settings\")\n", - "\n", - "try:\n", - " with exiftool.ExifTool(os.environ.get('exiftoolpath')) as exift:\n", - " print('Successfully executed exiftool.')\n", - "except Exception as e:\n", - " print(\"Exiftool isn't working. Double check that you've followed the instructions above.\")\n", - " print(\"The execption text below may help to find the source of the problem:\")\n", - " print()\n", - " print(e)\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Testing image reading and panel detection\n", - "\n", - "The above code checks for the proper libraries to be installed and verifies it can execute `exiftool`. This code opens an example image, reads the metadata, and then uses the `pyzbar` library to find a MicaSense panel in the image. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from micasense.image import Image\n", - "imagePath = os.path.join('.','data','0000SET','000')\n", - "imageName = glob.glob(os.path.join(imagePath,'IMG_0000_1.tif'))[0]\n", - "\n", - "img = Image(imageName)\n", - "img.plot_raw(figsize=(8.73,8.73));\n", - "\n", - "from micasense.panel import Panel\n", - "panel = Panel(img)\n", - "if not panel.panel_detected():\n", - " raise IOError(\"Panel Not Detected! Check your installation of pyzbar\")\n", - "else:\n", - " panel.plot(figsize=(8,8));\n", - "\n", - "print('Success! Now you are ready for Part 1 of the tutorial.')" - ] - }, - { - "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", - "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.9.2" - }, - "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": 4 -} diff --git a/MicaSense Image Processing Tutorial 1.ipynb b/MicaSense Image Processing Tutorial 1.ipynb deleted file mode 100755 index 1cc0c578..00000000 --- a/MicaSense Image Processing Tutorial 1.ipynb +++ /dev/null @@ -1,404 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# MicaSense RedEdge Image Processing Tutorial 1\n", - "\n", - "## Overview\n", - "\n", - "This tutorial assumes you have gone through the basic setup [here](./MicaSense Image Processing Setup.html) and your system is set up and ready to go.\n", - "\n", - "In this tutorial, we will walk through how to convert RedEdge data from raw images to radiace and then to reflectance. We will cover the tools required to do this, and walk through some of the basic image processing and radiometric conversions. \n", - "\n", - "### Opening an image with pyplot\n", - "\n", - "RedEdge 16-bit images can be read with pyplot directly into numpy arrays using the pyplot `imread` function or the matplotlib `imread` function, and then we can display the image inline using the `imshow` function of `matplotlib`.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import cv2\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import os,glob\n", - "import math\n", - "%matplotlib inline\n", - "\n", - "imagePath = os.path.join('.','data','0000SET','000')\n", - "imageName = os.path.join(imagePath,'IMG_0000_4.tif')\n", - "\n", - "# Read raw image DN values\n", - "# reads 16 bit tif - this will likely not work for 12 bit images\n", - "imageRaw=plt.imread(imageName)\n", - "\n", - "# Display the image\n", - "fig, ax = plt.subplots(figsize=(8,6))\n", - "ax.imshow(imageRaw, cmap='gray')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### MicaSense Utilities Module\n", - "\n", - "For many of the steps in the tutorial, we will use code from the MicaSense utilities module. The code is in the micasense directory and can be imported via normal python import commands using the syntax `import micasense` or `import micasense.submodule as short_name` for use in this and other scripts. While we will not cover all of the utility functions in this tutorial, they are available for reference and some will be used and discussed in future tutorials.\n", - "\n", - "### Adding a colorbar\n", - "\n", - "We will use start by using a plotting function in `micasense.plotutils` that adds a colorbar to the display, so that we can more easily see changes in the values in the images and also see the range of the image values after various conversions. This function also colorizes the grayscale images, so that changes can more easily be seen. Depending on your viewing style, you may prefer a different color map and you can also select that colormap here or browsing the colormaps on the [matplotlib site](https://matplotlib.org/users/colormaps.html). " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import micasense.plotutils as plotutils\n", - "\n", - "# Optional: pick a color map that fits your viewing style\n", - "# one of 'gray, viridis, plasma, inferno, magma, nipy_spectral'\n", - "plotutils.colormap('viridis'); \n", - "\n", - "fig = plotutils.plotwithcolorbar(imageRaw, title='Raw image values with colorbar')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reading RedEdge Metadata\n", - "\n", - "In order to perform various processing on the images, we need to read the metadata of each image. For this we use ExifTool. We can read standard image capture metadata such as location, UTC time, imager exposure and gain, but also RedEdge specific metadata which can make processing workflows easier. \n", - "\n", - "For example, each image contains a unique capture identifier. Capture identifiers are shared between all 5 images captured by RedEdge at the same moment, and can be used to unambiguously group images in post processing, regardless of how the images are named or stored on disk. Each image also contains a flight identifer which is the same for all images taken during a single power cycle of the camera. This can be used in post-processing workflows to group images and in many cases, more easily identify when the vehicle took off and landed." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import micasense.metadata as metadata\n", - "exiftoolPath = None\n", - "if os.name == 'nt':\n", - " exiftoolPath = os.environ.get('exiftoolpath')\n", - "# get image metadata\n", - "meta = metadata.Metadata(imageName, exiftoolPath=exiftoolPath)\n", - "cameraMake = meta.get_item('EXIF:Make')\n", - "cameraModel = meta.get_item('EXIF:Model')\n", - "firmwareVersion = meta.get_item('EXIF:Software')\n", - "bandName = meta.get_item('XMP:BandName')\n", - "print('{0} {1} firmware version: {2}'.format(cameraMake, \n", - " cameraModel, \n", - " firmwareVersion))\n", - "print('Exposure Time: {0} seconds'.format(meta.get_item('EXIF:ExposureTime')))\n", - "print('Imager Gain: {0}'.format(meta.get_item('EXIF:ISOSpeed')/100.0))\n", - "print('Size: {0}x{1} pixels'.format(meta.get_item('EXIF:ImageWidth'),meta.get_item('EXIF:ImageHeight')))\n", - "print('Band Name: {0}'.format(bandName))\n", - "print('Center Wavelength: {0} nm'.format(meta.get_item('XMP:CentralWavelength')))\n", - "print('Bandwidth: {0} nm'.format(meta.get_item('XMP:WavelengthFWHM')))\n", - "print('Capture ID: {0}'.format(meta.get_item('XMP:CaptureId')))\n", - "print('Flight ID: {0}'.format(meta.get_item('XMP:FlightId')))\n", - "print('Focal Length: {0}'.format(meta.get_item('XMP:FocalLength')))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Converting raw images to Radiance\n", - "\n", - "Ultimately most RedEdge users want to calibrate raw images from the camera into reflectance maps. This can be done using off-the-shelf software from third parties, but you are here because there is no fun in that! Along with this tutorial we have included some helper utilities that will handle much of this conversion for you, but here we will walk through a few of those functions to discuss what is happening inside.\n", - "\n", - "Any RedEdge workflow must include these common steps.\n", - "1. Un-bias images by accounting for the dark pixel offset\n", - "1. Compensate for imager-level effects\n", - "1. Compensate for optical chain effects\n", - "1. Normalize images by exposure and gain settings\n", - "1. Convert to a common unit system (radiance)\n", - "\n", - "All of these are handled by the `micasense.utils.raw_image_to_radiance(metadata, raw_image)` function. Let us take a look at that fuction in more detail.\n", - "\n", - "First, we get the darkPixel values. These values come from optically-covered pixels on the imager which are exposed at the same time as the image pixels. They measure the small amount of random charge generation in each pixel, independent of incoming light, which is common to all semiconductor imaging devices.\n", - "```python\n", - " blackLevel = np.array(meta.get_item('Exif.BlackLevel'))\n", - " darkLevel = blackLevel.mean()\n", - "```\n", - "Now, we get the imager-specific calibrations.\n", - "```python\n", - " a1, a2, a3 = meta.get_item('XMP:RadiometricCalibration')\n", - "```\n", - "We get the parameters of the optical chain (vignette) effects and create a vignette map. This map will be multiplied by the black-level corrected image values to reverse the darkening seen at the image corners. See the `vignette_map` function for the details of the vignette parameters and their use.\n", - "```python\n", - " V, x, y = vignette_map(meta, xDim, yDim)\n", - "```\n", - "Now we can calculate the imager-specfic radiometric correction function, which help to account for the radiometric inaccuracies of the CMOS imager pixels.\n", - "```python\n", - " # row gradient correction\n", - " R = 1.0 / (1.0 + a2 * y / exposureTime - a3 * y)\n", - "```\n", - "Finally, we apply these functions to the raw image to result in a corrected image\n", - "```python\n", - " # subtract the dark level and adjust for vignette and row gradient\n", - " L = V * R * (imageRaw - darkLevel)\n", - "```\n", - "Next, we get the exposure and gain settings (gain is represented in the photographic parameter ISO, with a base ISO of 100, so we divide the result to get a numeric gain).\n", - "```python\n", - " exposureTime = float(meta.get_item('EXIF:ExposureTime'))\n", - " gain = float(meta.get_item('EXIF:ISOSpeed'))/100.0 \n", - "```\n", - "Now that we have a corrected image, we can apply a conversion from calibrated digital number values to radiance units (W/m^2/nm/sr). Note that in this conversion, we need to normalize by the image bitdepth (2^16 for 16 bit images, 2^12 for 12-bit images), because the calibration coefficients are scaled to work with normalized input values.\n", - "```python\n", - " # apply the radiometric calibration - \n", - " # scale by the gain-exposure product and multiply with the radiometric calibration coefficient\n", - " bitsPerPixel = meta.get_item('EXIF:BitsPerSample')\n", - " dnMax = float(2**bitsPerPixel)\n", - " radianceImage = L.astype(float)/(gain * exposureTime)*a1/dnMax\n", - " ```\n", - " \n", - "For convenience, we have written the `raw_image_to_radiance` function to return the intermediate compensation images as well, so we can visualize them for the tutorial. These intermediate results are not required in most implementations and can be ommitted, if performance is a concern." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import micasense.utils as msutils\n", - "radianceImage, L, V, R = msutils.raw_image_to_radiance(meta, imageRaw)\n", - "plotutils.plotwithcolorbar(V,'Vignette Factor');\n", - "plotutils.plotwithcolorbar(R,'Row Gradient Factor');\n", - "plotutils.plotwithcolorbar(V*R,'Combined Corrections');\n", - "plotutils.plotwithcolorbar(L,'Vignette and row gradient corrected raw values');\n", - "plotutils.plotwithcolorbar(radianceImage,'All factors applied and scaled to radiance');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Convert radiance to reflectance\n", - "\n", - "Now that we have a flat and calibrated radiance image, we can convert into reflectance. To do this, we will use the radiance values of the panel image of known reflectance to determine a scale factor between radiance and reflectance.\n", - "\n", - "In this case, we have our MicaSense calibrated reflectance panel and it's known reflectance of 62% in the band of interest. We will extract the area of the image containing the lambertian panel, determine it's radiance to reflectance scale factor, and then scale the whole image by that factor to get a reflectance image." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "markedImg = radianceImage.copy()\n", - "ulx = 660 # upper left column (x coordinate) of panel area\n", - "uly = 490 # upper left row (y coordinate) of panel area\n", - "lrx = 840 # lower right column (x coordinate) of panel area\n", - "lry = 670 # lower right row (y coordinate) of panel area\n", - "cv2.rectangle(markedImg,(ulx,uly),(lrx,lry),(0,255,0),3)\n", - "\n", - "# Our panel calibration by band (from MicaSense for our specific panel)\n", - "panelCalibration = { \n", - " \"Blue\": 0.67, \n", - " \"Green\": 0.69, \n", - " \"Red\": 0.68, \n", - " \"Red edge\": 0.67, \n", - " \"NIR\": 0.61 \n", - "}\n", - "\n", - "# Select panel region from radiance image\n", - "panelRegion = radianceImage[uly:lry, ulx:lrx]\n", - "plotutils.plotwithcolorbar(markedImg, 'Panel region in radiance image')\n", - "meanRadiance = panelRegion.mean()\n", - "print('Mean Radiance in panel region: {:1.3f} W/m^2/nm/sr'.format(meanRadiance))\n", - "panelReflectance = panelCalibration[bandName]\n", - "radianceToReflectance = panelReflectance / meanRadiance\n", - "print('Radiance to reflectance conversion factor: {:1.3f}'.format(radianceToReflectance))\n", - "\n", - "reflectanceImage = radianceImage * radianceToReflectance\n", - "plotutils.plotwithcolorbar(reflectanceImage, 'Converted Reflectane Image');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In some cases we might notice that some parts of a converted reflectance image show reflectances above 1.0, or 100%, and wonder how this is possible. In fact, reflectances higher than 100% are normal in specific cases of [specular reflections](https://en.wikipedia.org/wiki/Specular_reflection). The panel area is a special material that reflects incident light equally well in all directions; however, some of the objects in a scene (especially man-made objects) instead reflect most incident light in one direction, more like a mirror. An example is the reflection of the sun off of the smooth surface of a car or the reflection off of a body of water.\n", - "\n", - "\n", - "Now we will extract the same region and verify the reflectance in that region is what we expect. In the process, we will blur and visualize the extracted area to look for any trends. The area should have a very consistent reflectance. If a gradient or a high standard deviation (>3% absolute reflectance) is noticed across the panel area it is possible that the panel was captured under inconsistent lighting conditions (e.g. next to a wall or vehicle) or it was captured too close to the edge of the image where the optical calibration is the least accurate." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "panelRegionRaw = imageRaw[uly:lry, ulx:lrx]\n", - "panelRegionRefl = reflectanceImage[uly:lry, ulx:lrx]\n", - "panelRegionReflBlur = cv2.GaussianBlur(panelRegionRefl,(55,55),5)\n", - "plotutils.plotwithcolorbar(panelRegionReflBlur, 'Smoothed panel region in reflectance image')\n", - "print('Min Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.min()))\n", - "print('Max Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.max()))\n", - "print('Mean Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.mean()))\n", - "print('Standard deviation in region: {:1.4f}'.format(panelRegionRefl.std()))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In this case the panel is less uniform than we would like, but we can also notice that the full color scale is only 4% of absolute reflectance. Likewise, it's well below the standard deviation threshold we have set of 3% absolute reflectance. This panel has seen over two years of hard field use, so it may be time for it to retire. \n", - "\n", - "Reasons for a high standard deviation across a panel can include panel contamination or inconsistent lighting across the panel due to environmental conditions. Based on the context of the image, it is also clear that the user is taking the panel image facing the sun, which can cast reflected light from the operator's clothing on the panel and contaminate results. For this reason it is always best to capture panel images in an open area and with the operator's back to the sun. \n", - "\n", - "### Undistorting images\n", - "\n", - "Finally, we need to remove lens distortion effects from images for some processing workflows, such as band-to-band image alignment. Generally for photogrammetry processes on raw (or radiance/reflectance) images, this step is not required, as the photogrammetry process will optimize a lens distortion model as part of it's bulk bundle adjustment. RedEdge has very low distortion lenses, so the changes to images in this step tend to be very small and noticeable only in pixels on the border of the image." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "# correct for lens distortions to make straight lines straight\n", - "undistortedReflectance = msutils.correct_lens_distortion(meta, reflectanceImage)\n", - "plotutils.plotwithcolorbar(undistortedReflectance, 'Undistorted reflectance image');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# In Practice\n", - "\n", - "Now that we can convert from raw RedEge images to reflectance, we will use these methods to convert an image taken during the same campaign to a reflectance image, and look at a few interesting areas of the image to validate our conversion." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flightImageName = os.path.join(imagePath,'IMG_0001_4.tif')\n", - "flightImageRaw=plt.imread(flightImageName)\n", - "plotutils.plotwithcolorbar(flightImageRaw, 'Raw Image')\n", - "\n", - "flightRadianceImage, _, _, _ = msutils.raw_image_to_radiance(meta, flightImageRaw)\n", - "flightReflectanceImage = flightRadianceImage * radianceToReflectance\n", - "flightUndistortedReflectance = msutils.correct_lens_distortion(meta, flightReflectanceImage)\n", - "plotutils.plotwithcolorbar(flightUndistortedReflectance, 'Reflectance converted and undistorted image');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Conclusion\n", - "\n", - "In this tutorial we have found that we can read MicaSense RedEdge images and their metadata, and use python and OpenCV to convert those images to radiance and then to reflectance using the standard scientific field method of imaging a lambertian reflector. We have corrected for both the electro-optical effects of the sensor and optical chain, as well as the incident light at the time of capture.\n", - "\n", - "In future tutorials, we will introduce the Downwelling Light Sensor (DLS) information into the calibration process in order to account for changing irradiance over time (e.g. such as clouds). However, since the panel method is straightforward and repeatable under constant illumination conditions, and is the standard scientific calibration method of surface reflectance, this process is useful and sufficient for many calibration needs.\n", - "\n", - "Looking for more? Try the second tutorial [here](./MicaSense%20Image%20Processing%20Tutorial%202.html)." - ] - }, - { - "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", - "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.1" - }, - "toc": { - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "toc_cell": false, - "toc_position": { - "height": "954px", - "left": "0px", - "right": "1123px", - "top": "107px", - "width": "249px" - }, - "toc_section_display": "block", - "toc_window_display": true - }, - "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 -} diff --git a/MicaSense Image Processing Tutorial 2.ipynb b/MicaSense Image Processing Tutorial 2.ipynb deleted file mode 100644 index 5c15319a..00000000 --- a/MicaSense Image Processing Tutorial 2.ipynb +++ /dev/null @@ -1,321 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Tutorial 2 - MicaSense library\n", - "\n", - "This tutorial assumes you have gone through the [basic setup](./Micasense Image Processing Setup.html) and builds on the basic radiance, irradiance, and reflectance concepts and code covered in the [first tutorial](./MicaSense Image Processing Tutorial 1.html). \n", - "\n", - "In this tutorial, we will cover usage of the MicaSense python library to access images and groups of images. Most of the processing details are hidden away in the library, but the library code is open and available in the git repository. \n", - "\n", - "# Library Components\n", - "\n", - "In the first tutorial, we introduced `micasense.utils` which provided some helper functions for single image manipulation, and `micasense.plotutils` which provided some plotting helpers.\n", - "\n", - "For this second tutorial, we are going to introduce the usage of the included micasense libraries for opening, converting, and displaying images. This will allow us to discuss and visualize results at a high level, while the underlying source code is available for those interested in the implementation details. In some cases, the libraries themselves may be enough to implement a custom workflow without the need to re-implement or translate the code to another system or language.\n", - "\n", - "The library code provides some basic classes to manage image data. At the highest level is the `ImageSet`, which is able to load a list of files or recursively search a whole directory into data structures which are easy to access and manipulate. `ImageSet`s are made up of `Capture`s, which hold the set of (usually 5) images as they are simultaneously gathered by the RedEdge camera. Within `Capture`s are `Image`s, which hold a single image file and allow easy access to the image metadata. The `Image` class also provides the ability to extract metadata from individual images and to convert individual images in similar ways to those described in the first tutorial. \n", - "\n", - "For the rest of this article, we will look at each of the objects available starting with the single `Image` object, and work our way up to the whole `ImageSet`. Each section in this article is standalone, and can be copied into another workbook or edited in place to explore more of the functions associated with that object." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## micasense.Image\n", - "\n", - "An image is the lowest level object. It represents the data in a single tiff file as taken by the camera. `Image` objects expose a set of data retrieval methods which provide access to raw, radiance, and reflectance corrected images, and to undistort any of those images. Note that when retrieving image data from an `Image` object, the data is stored internally in the object, increasing the object's memory footprint. If operating on a large number of images, it may be necessary to release this data memory after each image is processed to limit the program memory footprint. This can be done by calling the `Image.clear_image_data()` method." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import micasense.image as image\n", - "%matplotlib inline\n", - "\n", - "image_path = os.path.join('.','data','0000SET','000','IMG_0000_1.tif')\n", - "img = image.Image(image_path)\n", - "img.plot_raw();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Accessing `Image` Metadata\n", - "\n", - "Metadata for each image is available in the `Image.meta` parameter. This object is a `micasense.Metadata` object and can be accessed directly for image specific metadata extraction. Below, we print the same metadata values as we did in Tutorial #1, but using direct access to the `Metadata` object parameters.\n", - "\n", - "A notebook for experimenting with the `Image` class can be found [here](Images.html)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print('{0} {1} firmware version: {2}'.format(img.meta.camera_make(),\n", - " img.meta.camera_model(), \n", - " img.meta.firmware_version()))\n", - "print('Exposure Time: {0} seconds'.format(img.meta.exposure()))\n", - "print('Imager Gain: {0}'.format(img.meta.gain()))\n", - "print('Size: {0}x{1} pixels'.format(img.meta.image_size()[0],\n", - " img.meta.image_size()[1]))\n", - "print('Band Name: {0}'.format(img.meta.band_name()))\n", - "print('Center Wavelength: {0} nm'.format(img.meta.center_wavelength()))\n", - "print('Bandwidth: {0} nm'.format(img.meta.bandwidth()))\n", - "print('Capture ID: {0}'.format(img.meta.capture_id()))\n", - "print('Flight ID: {0}'.format(img.meta.flight_id()))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## micasense.Capture\n", - "\n", - "The `Capture` class is a container for `Image`s which allows access to metadata common to the group of images. The internal `Image` objects are accessible via the `capture.images` properties, and images in this list are kept sorted by the `band` property. Data which is different for each image can be accessed through composite methods, such as the `capture.dls_irradiance()` method, which returns a list of irradiances in band order." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os, glob\n", - "import micasense.capture as capture\n", - "\n", - "images_path = os.path.join('.','data','0000SET','000')\n", - "image_names = glob.glob(os.path.join(images_path,'IMG_0000_*.tif'))\n", - "cap = capture.Capture.from_filelist(image_names)\n", - "cap.plot_radiance();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Acessing `Capture` metadata\n", - "\n", - "Metadata which is common to all captures can be accessed via methods on the `Capture` object. Metadata which varies between the images of the capture, such as DLS information, is available as lists accessed from the capture object. \n", - "\n", - "