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", - "
\n", - " Note: The lists returned from metadata access on the `Capture` object are returned in `band_index` order. All images within a capture are sorted by the image `band_index`, and all lists adhere to this ordering. This ordering is consistent with the number at the end of each filename of a RedEdge image.\n", - "
\n", - "\n", - "Below we plot the raw and tilt compensated DLS irradiance by center wavelength and by band name.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "print(cap.band_names())\n", - "fig = plt.figure(figsize=(14,6))\n", - "plt.subplot(1,2,1)\n", - "plt.scatter(cap.center_wavelengths(), cap.dls_irradiance())\n", - "plt.ylabel('Irradiance $(W/m^2/nm)$')\n", - "plt.xlabel('Center Wavelength (nm)')\n", - "plt.subplot(1,2,2)\n", - "plt.scatter(cap.band_names(), [img.meta.exposure() for img in cap.images])\n", - "plt.xlabel('Band Names')\n", - "plt.ylim([0,2.5e-3])\n", - "plt.ylabel('Exposure Time (s)')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A notebook for experimenting with the `Capture` class can be found [here](Captures.html)." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## micasense.Panel\n", - "\n", - "The `Panel` class is a helper class which can automatically extract panel information from MicaSense calibrated reflectance panels by finding the QR code within an image and using the QR Code location and orientation information to find the lambertian panel area. The class then allows extraction of statistics from the panel area such as mean raw values, mean radiance, standard deviation, and the number of saturated pixels in the panel region. The panel object can be included standalone, or used within the context of a `Capture` object.\n", - "\n", - "
\n", - " Note: For the automatic panel QR code finding functions of the library to work, zbar and it's python bindings must be installed. We have made every effort to ensure this fails gracefully if zbar isn't available. Unfortunately zbar is only available using Python 2.7, not Python 3. If you're using Python 3.x, the code available in '/micasense/panel.py' shows how to find QR codes in images and to find the panel area from the QR location. We're currently looking for Python QR code finding options that work across platforms and Python versions, let us know if you have one that supports location!\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import os, glob\n", - "import micasense.image as image\n", - "import micasense.panel as panel\n", - "\n", - "image_path = os.path.join('.','data','0000SET','000','IMG_0000_1.tif')\n", - "img = image.Image(image_path)\n", - "# panelCorners - if we dont have zbar installed to scan the QR codes, detect panel manually and \n", - "panelCorners = [[[809,613],[648,615],[646,454],[808,452]],\n", - " [[772,623],[613,625],[610,464],[770,462]],\n", - " [[771,651],[611,653],[610,492],[770,490]],\n", - " [[829,658],[668,659],[668,496],[829,496]],\n", - " [[807,632],[648,634],[645,473],[805,471]]]\n", - "\n", - "pnl = panel.Panel(img,panelCorners = panelCorners[0])\n", - "print(\"Panel found: {}\".format(pnl.panel_detected()))\n", - "print(\"Panel serial: {}\".format(pnl.serial))\n", - "print(\"QR Code Corners:\\n{}\".format(pnl.qr_corners()))\n", - "mean, std, count, saturated_count = pnl.raw()\n", - "print(\"Panel mean raw pixel value: {}\".format(mean))\n", - "print(\"Panel raw pixel standard deviation: {}\".format(std))\n", - "print(\"Panel region pixel count: {}\".format(count))\n", - "print(\"Panel region saturated pixel count: {}\".format(count))\n", - "\n", - "pnl.plot();" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A notebook for experimenting with the `Panel` class can be found [here](Panels.html)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## micasense.ImageSet\n", - "\n", - "An `ImageSet` contains a group of `Capture`s. The captures can be loaded from image object, from a list of files, or by recursively searching a directory for images.\n", - "\n", - "Loading an `ImageSet` can be a time consuming process. It uses python multithreading under the hood to maximize cpu usage on multi-core machines." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import FloatProgress\n", - "from IPython.display import display\n", - "f = FloatProgress(min=0, max=1)\n", - "display(f)\n", - "def update_f(val):\n", - " f.value=val\n", - "\n", - "import micasense.imageset as imageset\n", - "import os\n", - "images_dir = os.path.join('.','data','0000SET')\n", - "\n", - "imgset = imageset.ImageSet.from_directory(images_dir, progress_callback=update_f)\n", - "\n", - "for cap in imgset.captures:\n", - " print (\"Opened Capture {} with bands {}\".format(cap.uuid,[str(band) for band in cap.band_names()]))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Extended ImageSet examples\n", - "\n", - "A large group of images captured 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, the [extended ImageSet example](./ImageSets.html) notebook provides more usages of ImageSet data.\n", - "\n", - "# Conclusion\n", - "\n", - "In this tutorial, we have introduced the MicaSense library and provided some examples of opening Images, Captures, and ImageSets, as well as detecting and extracting panel information from images.\n", - "\n", - "The next tutorial covers basic usage of DLS information, and is available [here](./MicaSense%20Image%20Processing%20Tutorial%203.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": 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 3.ipynb b/MicaSense Image Processing Tutorial 3.ipynb deleted file mode 100644 index 74718c30..00000000 --- a/MicaSense Image Processing Tutorial 3.ipynb +++ /dev/null @@ -1,272 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Tutorial 3 - DLS Sensor Basic Usage\n", - "\n", - "\n", - "This tutorial will walk through how to convert RedEdge data from raw images to radiace and then use the DLS information in the image metadata to convert that to reflectance.\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Update for DLS2\n", - "\n", - "This tutorial was originally written as a guide to basic processing of DLS1 data for RedEdge-3 and RedeEdge-M. Since October 2018, Altum and RedEdge-MX have been shipped standard with DLS2. DLS1 is a small, red square with a single diffuser on top. DLS2 is a larger, thinner black rectangle with two small diffusers on top and 8 small diffusers on different up-looking faces.\n", - "\n", - "For DLS2 data, we recommend using the `Capture` class to access the `Capture.dls_irradiance()` method. This will provide a compensated horizontal irradiance useful for radiometrically correcting imagery. We leave the below intact for legacy users and as a tutorial on remote sensing calibrations and terms. \n", - "\n", - "For RedEdge and Altum with up-to-date firmware, the `EXIF:HorizontalIrradiance` tag and the appropriate scale factor will provide the necessary information for reflectance compensation. If necessary, this can be further field-calibrated by using a pane image to tie the DLS2 and camera calibrations together in-situ." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Remote Sensing Theory\n", - "\n", - "In [Tutorial 1](./MicaSense%20Image%20Processing%20Tutorial%201.html) we covered conversion of the image radiance to reflecance using the average radiance from an area of a specially made lambertian reflectance panel. As we did not get very far into remote sensing theory, we will cover that some here as we move on to the more complex problem of using the data from the RedEdge Downwelling Light Sensor (DLS). To get started, we will review some of the fundamental measurements in remote sensing.\n", - "\n", - "## Definition of Terms\n", - "\n", - "\n", - "The DLS measures the energy coming from the sky, while the camera measures the energy coming from the fields below. The property that we are generally looking to estimate in agricultural remote sensing is the surface reflectance of the plants and ground below the sensor. \n", - "\n", - "The MicaSense RedEdge measures spectral radiance, a measure of the power incident on the CMOS image sensor at each pixel, and has units power per unit area per unit bandwidth per unit solid angle (W/m^2/nm/sr). We will call the radiance image L.\n", - "\n", - "The MicaSense Downwelling Light Sensor (DLS) measures incoming spectral irradiance on the top surface of the device for the same 5 spectral bands as the camera. Irradiance has the units of power per unit area per unit bandwidth, or Watts per meter squared per nanometer (W/m^2/nm). We will call this downwelling quantity I.\n", - "\n", - "If we look closely at the units, they differ by the last term - the unit over solid angle. This is because each pixel of the camera sensor only sees a specific solid angle, while the dls detector has a specific area. \n", - "\n", - "\n", - "## Understanding the measurements\n", - "\n", - "In order to use downwelling irradiance and upwelling radiance to determine reflectance, we simply need to divide pi times the irradiance by the radiance. Our units cancel out, and we are left with the unitless quantity of reflectance. However, there is one caveat. In most systems, the DLS is not held level to the ground at all times, so the amount of light that shines on it's surface changes as the aircraft moves about. In oder for that simple math to work, it's necessary that the irradiance and radiance be compensated for a very important property: the illumination angle. \n", - "\n", - "Both the radiance and irradiance are sensitive to changes in the direction of the lighting. Both quantites decrease with the cosine of the angle from the surface normal, eventually reaching zero when the light source is parallel to the surface. So when measuring irradiance, this effect needs to be corrected in order to estimate the actual irradiance.\n", - "\n", - "## Different types of irradiance\n", - "\n", - "To complicate matters further, the light reaching the DLS and the ground is not coming from a single source, even under sunny conditions. The direct sunlight provides what we will call the direct source, resulting in a direct irradiance `I_direct`. The rest of the (usually blue) sky provides what we term the indirect, or diffuse, irradiance which we will call `I_diffuse`. Under even the sunniest conditions with the summer sun directly overhead, only about 85% of the light is direct, while the other 15% is diffuse, which we will express as the ratio 6:1. Under thin cloud conditions, this ratio drops to close 1:1, meaning half of the measured irradiance is coming from parts of the sky other than the direction of the sun. The ratio goes even lower under fully cloudy skies, becoming fully diffuse.\n", - "\n", - "As we ultimately want to estimate the time-varying irradiance, we will need to estimate both the direct and diffuse irradiance terms as well as correct for the cosine effect of the DLS sensor itself.\n", - "\n", - "## Sensor-specific challenges\n", - "\n", - "There's one final piece of theory required that pertains to irradiance sensing, and that is the Fresnel Effect. The Fresnel Effect describes a physical property of the diffuser that gathers light on top of the DLS. This effect results from the fact that much of the light striking the diffuser is reflected, and the amount reflected changes with the angle of the incident light. So, instead of our diffuser behaving as a perfect cosine detector, we need to also compensate for the material imperfections.\n", - "\n", - "## DLS pose estimation\n", - "\n", - "An onboard orientation sensor provides a fairly rough version of the earth-fixed orientation of the DLS. The quality of this orientation varies based on the quality of the sensor installation and calibration, but in even the best cases it has a few degrees of error. If the DLS and the camera are mounted fixed to one another, the photogrammetric workflow can be used to augment the orientation of the irradiance sensor. Our examples here will use the sensor measurements directly, but significant improvements can be made using the very precise angles photorammetry can provide.\n", - "\n", - "## Radiometric calibration options\n", - "\n", - "As we've discussed in the first tutorial, the method of using a calibrated reflectance panel is the scientific standard for remote sensing measurements. Spectroradiometers which cost tens of thousands of dollars are calibrated in this method in the field. Panels images may also be used both before and after each flight to track slowly changing solar irradiance. In many cases of direct sunlight with no clouds, panel images can provide the best results when taken before and after each flight. As the diffuse lighting increases, the utility of the DLS increases as well.\n", - "\n", - "In cases of more varied irradiance, the DLS can be used for irradiance estimation, using only the internal pose estimation. The DLS irradiance can be be used in conjunction with panels to ensure that the DLS and panel measurements match at the time of panel capture. In order for this to be effective, it's important that the DLS and panel are exposed to the same irradiance. This is best accomplished by ensuring panel images are taken with the operator's back to the sun and the operator shadow next to the panel. The aircraft and operator should block as little of the sky as possible.\n", - "\n", - "Finally, in fixed mount systems, the photogrammetric angles can be used, along with the offset between the camera and DLS sensors, to most accurately estimate irradiance. This can again be coupled with one or more panel measurements for best radiometric results." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Putting it all Together\n", - "\n", - "Below, we will use the `micasense.Capture` class to load a set of images and use tilt-compensated DLS irradiance data to roughly translate those images to reflectance. We will then compare these to the known reflectances of a lambertian panel, and develop a site-specific compensation factor to tie together the DLS irradiance and camera radiance measurements through the panel.\n", - "\n", - "## Accessing and processing DLS data\n", - "\n", - "Below, we show how to access DLS data and compute the minimum required corrections to the DLS data for use. The code in this cell is present in the `micasense.dls` library, and is used by the `Image` and `Capture` class to expose corrected DLS data. If your intent is to simply use this code to correct your data, you can skip this section and access pose-corrected DLS data through the `capture.dls_irradiance()` method.\n", - "\n", - "For the purpose of this tutorial, we assume a constant irradiance (over time) and a ratio of direct to diffuse irradiance of 6:1. These assumptions only hold approximately, and only for clear sky. Cloud cover can lead to dramatically different results. \n", - "\n", - "## Loading Captures" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import micasense.dls as dls\n", - "\n", - "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", - "# set panel corners manually if zbar is not installed\n", - "panelCorners = [[[809,613],[648,615],[646,454],[808,452]],\n", - " [[772,623],[613,625],[610,464],[770,462]],\n", - " [[771,651],[611,653],[610,492],[770,490]],\n", - " [[829,658],[668,659],[668,496],[829,496]],\n", - " [[807,632],[648,634],[645,473],[805,471]]]\n", - "\n", - "cap.set_panelCorners(panelCorners)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Computing Solar orientation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Define DLS sensor orientation vector relative to dls pose frame\n", - "dls_orientation_vector = np.array([0,0,-1])\n", - "# compute sun orientation and sun-sensor angles\n", - "(\n", - " sun_vector_ned, # Solar vector in North-East-Down coordinates\n", - " sensor_vector_ned, # DLS vector in North-East-Down coordinates\n", - " sun_sensor_angle, # Angle between DLS vector and sun vector\n", - " solar_elevation, # Elevation of the sun above the horizon\n", - " solar_azimuth, # Azimuth (heading) of the sun\n", - ") = dls.compute_sun_angle(cap.location(),\n", - " cap.dls_pose(),\n", - " cap.utc_time(),\n", - " dls_orientation_vector)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Correcting DLS readings for orientations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Since the diffuser reflects more light at shallow angles than at steep angles,\n", - "# we compute a correction for this\n", - "fresnel_correction = dls.fresnel(sun_sensor_angle)\n", - "\n", - "# Now we can correct the raw DLS readings and compute the irradiance on level ground\n", - "dls_irradiances = []\n", - "center_wavelengths = []\n", - "for img in cap.images:\n", - " dir_dif_ratio = 6.0\n", - " percent_diffuse = 1.0/dir_dif_ratio\n", - " # measured Irradiance / fresnelCorrection\n", - " sensor_irradiance = img.spectral_irradiance / fresnel_correction\n", - " untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(sun_sensor_angle))\n", - " # compute irradiance on the ground using the solar altitude angle\n", - " dls_irr = untilted_direct_irr * (percent_diffuse + np.sin(solar_elevation))\n", - " dls_irradiances.append(dls_irr)\n", - " center_wavelengths.append(img.center_wavelength)\n", - "\n", - "import matplotlib.pyplot as plt\n", - "plt.scatter(center_wavelengths,dls_irradiances)\n", - "plt.xlabel('Wavelength (nm)')\n", - "plt.ylabel('Irradiance ($W/m^2/nm$)')\n", - "plt.show();\n", - "\n", - "cap.plot_undistorted_reflectance(dls_irradiances)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Tying the camera and DLS together via Panels\n", - "\n", - "Finally, if we want to use the panel as the ultimate reference tying the panel and the DLS readings together with a known lambertian reflector, we can calculate a DLS factor which we can then apply to DLS derived irradiances over the whole flight. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import math\n", - "\n", - "panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] #RedEdge band_index order\n", - "\n", - "panel_radiances = np.array(cap.panel_radiance())\n", - "irr_from_panel = math.pi * panel_radiances / panel_reflectance_by_band\n", - "dls_correction = irr_from_panel/dls_irradiances\n", - "cap.plot_undistorted_reflectance(dls_irradiances*dls_correction)\n", - "\n", - "plt.scatter(cap.center_wavelengths(), cap.panel_reflectance())\n", - "plt.title(\"Panel Reflectances\")\n", - "plt.xlabel(\"Wavelength (nm)\")\n", - "plt.ylabel(\"Reflectance\")\n", - "plt.show()\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "# Future improvements\n", - "\n", - "While the need to compensate the DLS for movement increases the complexity of the problem, it also provides us more information than might be available for other diffuser types. This is because as the DLS moves about on the aircraft, it provides many measurements per second that should all have consistent effects assuming the light isn't changing over those very fast measurements. This movement can then be used to develop an irradiance series which uses all of the available information from the on board sensors. A future tutorial will provide an overview of some methods for improving irradiance sensing further." - ] - }, - { - "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": "1002px", - "left": "0px", - "right": "1117px", - "top": "107px", - "width": "212px" - }, - "toc_section_display": "block", - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/Panels.ipynb b/Panels.ipynb deleted file mode 100644 index 6568f554..00000000 --- a/Panels.ipynb +++ /dev/null @@ -1,151 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Panels\n", - "This notebook shows usage for the Panel class. This type is useful for detecting MicaSense calibration panels and extracting information about the lambertian panel surface." - ] - }, - { - "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", - "from micasense.image import Image\n", - "from micasense.panel import Panel\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", - "panel = Panel(img)\n", - "\n", - "if not panel.panel_detected():\n", - " raise IOError(\"Panel Not Detected!\")\n", - " \n", - "print(\"Detected panel serial: {}\".format(panel.serial))\n", - "mean, std, num, sat_count = panel.raw()\n", - "print(\"Extracted Panel Statistics:\")\n", - "print(\"Mean: {}\".format(mean))\n", - "print(\"Standard Deviation: {}\".format(std))\n", - "print(\"Panel Pixel Count: {}\".format(num))\n", - "print(\"Saturated Pixel Count: {}\".format(sat_count))\n", - "\n", - "panel.plot();" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "imagePath = os.path.join('.','data','ALTUM1SET','000')\n", - "imageName = glob.glob(os.path.join(imagePath,'IMG_0000_1.tif'))[0]\n", - "\n", - "img = Image(imageName)\n", - "if img.auto_calibration_image:\n", - " print(\"Found automatic calibration image\")\n", - "panel = Panel(img)\n", - "\n", - "if not panel.panel_detected():\n", - " raise IOError(\"Panel Not Detected!\")\n", - " \n", - "print(\"Detected panel serial: {}\".format(panel.serial))\n", - "mean, std, num, sat_count = panel.raw()\n", - "print(\"Extracted Panel Statistics:\")\n", - "print(\"Mean: {}\".format(mean))\n", - "print(\"Standard Deviation: {}\".format(std))\n", - "print(\"Panel Pixel Count: {}\".format(num))\n", - "print(\"Saturated Pixel Count: {}\".format(sat_count))\n", - "\n", - "panel.plot();" - ] - }, - { - "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.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/data/0000SET/000/IMG_0000_1.tif b/data/0000SET/000/IMG_0000_1.tif deleted file mode 100644 index c64d7a57..00000000 --- a/data/0000SET/000/IMG_0000_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:248f42023cda80cae83944027b96d284d9b5f4573892dc597e99677ffe472eef -size 2465274 diff --git a/data/0000SET/000/IMG_0000_2.tif b/data/0000SET/000/IMG_0000_2.tif deleted file mode 100644 index a347f9a1..00000000 --- a/data/0000SET/000/IMG_0000_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e521e805839c129ea6b23c98e785a9f63486db042b05e6fdc9c248723909ecea -size 2465278 diff --git a/data/0000SET/000/IMG_0000_3.tif b/data/0000SET/000/IMG_0000_3.tif deleted file mode 100644 index e7625663..00000000 --- a/data/0000SET/000/IMG_0000_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:110d36022643142501a04f8ceac1309abf43bd03fe39cde9ee0101676986e0eb -size 2465274 diff --git a/data/0000SET/000/IMG_0000_4.tif b/data/0000SET/000/IMG_0000_4.tif deleted file mode 100644 index 8d691162..00000000 --- a/data/0000SET/000/IMG_0000_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:488fcc857be9583e69f533c010a82eee7e96dc690e266678f896f78f8d4e84ea -size 2465274 diff --git a/data/0000SET/000/IMG_0000_5.tif b/data/0000SET/000/IMG_0000_5.tif deleted file mode 100644 index 1d2d4c56..00000000 --- a/data/0000SET/000/IMG_0000_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:47d7a1a0c2eb289ce8ee1069c861e2a9838c7441d5038d6e4c9c9e37af1299fc -size 2465280 diff --git a/data/0000SET/000/IMG_0001_1.tif b/data/0000SET/000/IMG_0001_1.tif deleted file mode 100644 index 56d43656..00000000 --- a/data/0000SET/000/IMG_0001_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ee210953eefc85e3551b143f91b65c058743afdf4a2ff5a8cac98350cb8b7c8e -size 2465280 diff --git a/data/0000SET/000/IMG_0001_2.tif b/data/0000SET/000/IMG_0001_2.tif deleted file mode 100644 index 19d2f54b..00000000 --- a/data/0000SET/000/IMG_0001_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:d1b25a7317b45ae5d64c8111048faac74ed77ddb2b29533b3b50218a6c1eccac -size 2465282 diff --git a/data/0000SET/000/IMG_0001_3.tif b/data/0000SET/000/IMG_0001_3.tif deleted file mode 100644 index d0564132..00000000 --- a/data/0000SET/000/IMG_0001_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:048230b498e5f8d8e68935e4ecc64de99a3899b7900e721dcd876034fdd0ad93 -size 2465276 diff --git a/data/0000SET/000/IMG_0001_4.tif b/data/0000SET/000/IMG_0001_4.tif deleted file mode 100644 index 4ee5e1aa..00000000 --- a/data/0000SET/000/IMG_0001_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3581de2646fcf28b42a059551a1938ddda3ba744f954d6246748a4af85e8ecf2 -size 2465278 diff --git a/data/0000SET/000/IMG_0001_5.tif b/data/0000SET/000/IMG_0001_5.tif deleted file mode 100644 index ebd96c48..00000000 --- a/data/0000SET/000/IMG_0001_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6240d3c9fd8e89ee57c6d06ff837a38b40e12c7a3456dc19ac6bbcb5bf4bd833 -size 2465282 diff --git a/data/0001SET/000/IMG_0002_4.tif b/data/0001SET/000/IMG_0002_4.tif deleted file mode 100755 index 732f64f4..00000000 --- a/data/0001SET/000/IMG_0002_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1187b500f0196ce91e3b88a110e154b848a68d2a29bae2d9f69a30bc968bc88e -size 2465342 diff --git a/data/0001SET/000/IMG_0003_1.tif b/data/0001SET/000/IMG_0003_1.tif deleted file mode 100755 index 83c26d36..00000000 --- a/data/0001SET/000/IMG_0003_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ad25ef1f833b633a6652ed798d324efa1b7220569f71dc2419bcfe23108c2a4b -size 2464016 diff --git a/data/0002SET/000/IMG_0000_1.tif b/data/0002SET/000/IMG_0000_1.tif deleted file mode 100644 index b2128f10..00000000 --- a/data/0002SET/000/IMG_0000_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e0ef38734a94aae45347f96512d1d084d470956e805c0907a8d7aa659d289ce6 -size 2466318 diff --git a/data/0002SET/000/IMG_0000_2.tif b/data/0002SET/000/IMG_0000_2.tif deleted file mode 100644 index 421c52e0..00000000 --- a/data/0002SET/000/IMG_0000_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:aa00bc26b2037d67f1de7c4618b9d3ac6fc1bea5e22851dcf13a3cdcdad8a400 -size 2466308 diff --git a/data/0002SET/000/IMG_0000_3.tif b/data/0002SET/000/IMG_0000_3.tif deleted file mode 100644 index b4915698..00000000 --- a/data/0002SET/000/IMG_0000_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e41e98bea9994eff7959f4e93c134f0185f8287e2cda1bd19493528de3610b6f -size 2466300 diff --git a/data/0002SET/000/IMG_0000_4.tif b/data/0002SET/000/IMG_0000_4.tif deleted file mode 100644 index 0e2702e2..00000000 --- a/data/0002SET/000/IMG_0000_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:008b19f85db896219d662eb88434103e8393af09e43fbbe00f0be0fa40e848b8 -size 2466336 diff --git a/data/0002SET/000/IMG_0000_5.tif b/data/0002SET/000/IMG_0000_5.tif deleted file mode 100644 index df248a6b..00000000 --- a/data/0002SET/000/IMG_0000_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1c149e3bc4b0edc16c351c46922aec0aa42c14a820e00de818edf7c9a6252169 -size 2466332 diff --git a/data/10BANDSET/000/IMG_0000_1.tif b/data/10BANDSET/000/IMG_0000_1.tif deleted file mode 100644 index 6aed3457..00000000 --- a/data/10BANDSET/000/IMG_0000_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:a40134121e1b10440279b37da6f489a70046ee8bb1b570ab32e486882f2b3848 -size 2466196 diff --git a/data/10BANDSET/000/IMG_0000_10.tif b/data/10BANDSET/000/IMG_0000_10.tif deleted file mode 100644 index 7f6a9b04..00000000 --- a/data/10BANDSET/000/IMG_0000_10.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:600a305204cbab2b7fe0dc5c6398ee68e8c7b77605aa4de3d300124ff09d30de -size 2465694 diff --git a/data/10BANDSET/000/IMG_0000_2.tif b/data/10BANDSET/000/IMG_0000_2.tif deleted file mode 100644 index 8a6c1347..00000000 --- a/data/10BANDSET/000/IMG_0000_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:72991a75c5e0804744cb16b97ec88da24391602203e19dd0a05140b059ef2694 -size 2466188 diff --git a/data/10BANDSET/000/IMG_0000_3.tif b/data/10BANDSET/000/IMG_0000_3.tif deleted file mode 100644 index 5026afcc..00000000 --- a/data/10BANDSET/000/IMG_0000_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:2c5bfcf7596d76a853ff5e945fa83b98713f03d1e4ec61f0a693c42c81cf9a17 -size 2466188 diff --git a/data/10BANDSET/000/IMG_0000_4.tif b/data/10BANDSET/000/IMG_0000_4.tif deleted file mode 100644 index ce9b1242..00000000 --- a/data/10BANDSET/000/IMG_0000_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6dd4f58759ef4b2627728eb4c28adca711ece017223f8cca64eb6bfd2940cc5f -size 2466196 diff --git a/data/10BANDSET/000/IMG_0000_5.tif b/data/10BANDSET/000/IMG_0000_5.tif deleted file mode 100644 index e2ea0e58..00000000 --- a/data/10BANDSET/000/IMG_0000_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8f6d0e3802c78a4ce5e107f26c97cc4ad1e5fe6a87c78955a4fb41e19382f735 -size 2466198 diff --git a/data/10BANDSET/000/IMG_0000_6.tif b/data/10BANDSET/000/IMG_0000_6.tif deleted file mode 100644 index 0ddb9965..00000000 --- a/data/10BANDSET/000/IMG_0000_6.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1e37f9786c8a51f313a6437e6ddd3fcee6f198bedeaf53ffb53cb4210ca6b1b0 -size 2465698 diff --git a/data/10BANDSET/000/IMG_0000_7.tif b/data/10BANDSET/000/IMG_0000_7.tif deleted file mode 100644 index 6cba9be3..00000000 --- a/data/10BANDSET/000/IMG_0000_7.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:493bdd84bba111dd35ddc094d620c3acc01af24907a65191b0625c237352accf -size 2465714 diff --git a/data/10BANDSET/000/IMG_0000_8.tif b/data/10BANDSET/000/IMG_0000_8.tif deleted file mode 100644 index 52440229..00000000 --- a/data/10BANDSET/000/IMG_0000_8.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:69387338b89dfab6723f0b5d017087fa566697563a7691d51527a80f93165935 -size 2465648 diff --git a/data/10BANDSET/000/IMG_0000_9.tif b/data/10BANDSET/000/IMG_0000_9.tif deleted file mode 100644 index f45e9a6d..00000000 --- a/data/10BANDSET/000/IMG_0000_9.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:500d7d89c4621373b6c7ca8a4779b4171f505da41543c301d49cf1d5a50157cd -size 2465674 diff --git a/data/10BANDSET/000/IMG_0431_1.tif b/data/10BANDSET/000/IMG_0431_1.tif deleted file mode 100644 index a53a4b3a..00000000 --- a/data/10BANDSET/000/IMG_0431_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:20ae60bea30d14fcaa51b6ed6e729f5e4e47ecb9540e22d5c9ace3d3faa96074 -size 2465850 diff --git a/data/10BANDSET/000/IMG_0431_10.tif b/data/10BANDSET/000/IMG_0431_10.tif deleted file mode 100644 index 6555f46a..00000000 --- a/data/10BANDSET/000/IMG_0431_10.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1a21d4c66d8c0630bac9d959298cd3f9e29393d57c80a82715c52eee43247e3c -size 2465332 diff --git a/data/10BANDSET/000/IMG_0431_2.tif b/data/10BANDSET/000/IMG_0431_2.tif deleted file mode 100644 index 7ec56290..00000000 --- a/data/10BANDSET/000/IMG_0431_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ac3b8d2d9a45f27e18c85db103cc95558e1deb03be813d25b3a210700c30a506 -size 2465830 diff --git a/data/10BANDSET/000/IMG_0431_3.tif b/data/10BANDSET/000/IMG_0431_3.tif deleted file mode 100644 index c80d0651..00000000 --- a/data/10BANDSET/000/IMG_0431_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:74008850326fc703f54ca46b227e8f106c461e966c9f8e963380cac3cc667ae2 -size 2465832 diff --git a/data/10BANDSET/000/IMG_0431_4.tif b/data/10BANDSET/000/IMG_0431_4.tif deleted file mode 100644 index 7c930105..00000000 --- a/data/10BANDSET/000/IMG_0431_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:56dbbb7dddb15c10f2ed8beef4e5896be69bfd726d2f235ab5ab26485ed1c0e8 -size 2465840 diff --git a/data/10BANDSET/000/IMG_0431_5.tif b/data/10BANDSET/000/IMG_0431_5.tif deleted file mode 100644 index 9ce7c07f..00000000 --- a/data/10BANDSET/000/IMG_0431_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8cd66d530bdf11b68cc41d1caddfdba8abf4c4ff31f36da9ae5971f335e17989 -size 2465844 diff --git a/data/10BANDSET/000/IMG_0431_6.tif b/data/10BANDSET/000/IMG_0431_6.tif deleted file mode 100644 index 2d5bc074..00000000 --- a/data/10BANDSET/000/IMG_0431_6.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3d2e89fbea0515aca8c49302e889064981baa418f38f9910e9ba184707bb2e01 -size 2465338 diff --git a/data/10BANDSET/000/IMG_0431_7.tif b/data/10BANDSET/000/IMG_0431_7.tif deleted file mode 100644 index 73ca853d..00000000 --- a/data/10BANDSET/000/IMG_0431_7.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ba4699bfd1dd07946b7d87da81b6fb0bcfa9dd03e1233e5f47e8e3c9890de26c -size 2465350 diff --git a/data/10BANDSET/000/IMG_0431_8.tif b/data/10BANDSET/000/IMG_0431_8.tif deleted file mode 100644 index 73fe1fe0..00000000 --- a/data/10BANDSET/000/IMG_0431_8.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:24f6f8657e2ef4b4589764044eecc608992bfe5950f984bd7cb27c1346678382 -size 2465286 diff --git a/data/10BANDSET/000/IMG_0431_9.tif b/data/10BANDSET/000/IMG_0431_9.tif deleted file mode 100644 index d660227b..00000000 --- a/data/10BANDSET/000/IMG_0431_9.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:0fb1196332670b30083c6bf468fa79e4fb15eaf0f7af5c24c7d9db9656b8818c -size 2465314 diff --git a/data/ALTUM0SET/000/IMG_0000_1.tif b/data/ALTUM0SET/000/IMG_0000_1.tif deleted file mode 100644 index 8511c057..00000000 --- a/data/ALTUM0SET/000/IMG_0000_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3410bf21c9fd553937f5e0cd87e57fe7c145b1b7e3dbe8ccf4d4ee0063cc2d3b -size 6381820 diff --git a/data/ALTUM1SET/000/IMG_0000_1.tif b/data/ALTUM1SET/000/IMG_0000_1.tif deleted file mode 100644 index 77dbd14c..00000000 --- a/data/ALTUM1SET/000/IMG_0000_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:9a2abd9511ec95130dc1550dae07f49fa350a26db98d8596ec245a2a5a036e17 -size 6382362 diff --git a/data/ALTUM1SET/000/IMG_0000_2.tif b/data/ALTUM1SET/000/IMG_0000_2.tif deleted file mode 100644 index af4d833c..00000000 --- a/data/ALTUM1SET/000/IMG_0000_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:02ae3e49b8a7c04537e6e8a73e574cbe0f1a5d2ac429b7afb06f4ffbfebaf6d9 -size 6382386 diff --git a/data/ALTUM1SET/000/IMG_0000_3.tif b/data/ALTUM1SET/000/IMG_0000_3.tif deleted file mode 100644 index c130c9a2..00000000 --- a/data/ALTUM1SET/000/IMG_0000_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:bf09f520f8c79c5e4dd946cb62c4ae75ef729cb4874d60392aac09923734a2a4 -size 6382404 diff --git a/data/ALTUM1SET/000/IMG_0000_4.tif b/data/ALTUM1SET/000/IMG_0000_4.tif deleted file mode 100644 index 1d1fde24..00000000 --- a/data/ALTUM1SET/000/IMG_0000_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:deb3e1f4af5a03b8335db264cfb86640605b0cff33847545463f57a480a5a5fc -size 6382362 diff --git a/data/ALTUM1SET/000/IMG_0000_5.tif b/data/ALTUM1SET/000/IMG_0000_5.tif deleted file mode 100644 index 1dde1d6d..00000000 --- a/data/ALTUM1SET/000/IMG_0000_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:768cf0987000cdcd7eee684f659ef9e230fe817f14dc682ff18a82c093786315 -size 6382376 diff --git a/data/ALTUM1SET/000/IMG_0000_6.tif b/data/ALTUM1SET/000/IMG_0000_6.tif deleted file mode 100644 index 7876a307..00000000 --- a/data/ALTUM1SET/000/IMG_0000_6.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:7ac6cec517accf2618fc57b9096a39d78fa7810b805ea6e75ce05ad3d6fcdeae -size 44730 diff --git a/data/ALTUM1SET/000/IMG_0008_1.tif b/data/ALTUM1SET/000/IMG_0008_1.tif deleted file mode 100644 index 0424f862..00000000 --- a/data/ALTUM1SET/000/IMG_0008_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3df18c9843f0dd4fcf5174a876988d72a977bc3dd9e14f7c69aa28e158fc18f8 -size 6382004 diff --git a/data/ALTUM1SET/000/IMG_0008_2.tif b/data/ALTUM1SET/000/IMG_0008_2.tif deleted file mode 100644 index 109c72ea..00000000 --- a/data/ALTUM1SET/000/IMG_0008_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:29a268948de0a4ca87d27d5ca4e98a8abf6dc8a58ced04f438bd5b9805f98a30 -size 6382002 diff --git a/data/ALTUM1SET/000/IMG_0008_3.tif b/data/ALTUM1SET/000/IMG_0008_3.tif deleted file mode 100644 index 79c2a06c..00000000 --- a/data/ALTUM1SET/000/IMG_0008_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:a7da92c38d7722a06c6d301273407ce2a08a72d4d95e0f1f6372ba732dbb3f8f -size 6382034 diff --git a/data/ALTUM1SET/000/IMG_0008_4.tif b/data/ALTUM1SET/000/IMG_0008_4.tif deleted file mode 100644 index 9a66e91d..00000000 --- a/data/ALTUM1SET/000/IMG_0008_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8bb77e029f50a7600f4ff719ab6f8b034a55e19dca617780ddf753d323f16ac8 -size 6381992 diff --git a/data/ALTUM1SET/000/IMG_0008_5.tif b/data/ALTUM1SET/000/IMG_0008_5.tif deleted file mode 100644 index 8714b94c..00000000 --- a/data/ALTUM1SET/000/IMG_0008_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ff8f028694602573268c97e025203596c5feca5435f982826bba1abf1e08b16c -size 6382006 diff --git a/data/ALTUM1SET/000/IMG_0008_6.tif b/data/ALTUM1SET/000/IMG_0008_6.tif deleted file mode 100644 index a2ee8687..00000000 --- a/data/ALTUM1SET/000/IMG_0008_6.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:cf26fac7549bd740b68ad43c19eb8a65be0360e4391e6762635cbbec7752266b -size 44730 diff --git a/data/ALTUM1SET/000/IMG_0245_1.tif b/data/ALTUM1SET/000/IMG_0245_1.tif deleted file mode 100644 index 3697e4e8..00000000 --- a/data/ALTUM1SET/000/IMG_0245_1.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:35aedd3a3700eedf6ffd5b6c4cd0c6469d4c1ca4801041ca79d9765d87123616 -size 6381990 diff --git a/data/ALTUM1SET/000/IMG_0245_2.tif b/data/ALTUM1SET/000/IMG_0245_2.tif deleted file mode 100644 index 23a12050..00000000 --- a/data/ALTUM1SET/000/IMG_0245_2.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:51cab65288a589b383050a829fba01e027ae66bb3aecbc6926b17882a8601c37 -size 6382002 diff --git a/data/ALTUM1SET/000/IMG_0245_3.tif b/data/ALTUM1SET/000/IMG_0245_3.tif deleted file mode 100644 index 89a8e2e5..00000000 --- a/data/ALTUM1SET/000/IMG_0245_3.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:d34872e318624b168882b3180b69f16bb0b625fbf89176f0a49d1255a29de21d -size 6382032 diff --git a/data/ALTUM1SET/000/IMG_0245_4.tif b/data/ALTUM1SET/000/IMG_0245_4.tif deleted file mode 100644 index 51a46e6d..00000000 --- a/data/ALTUM1SET/000/IMG_0245_4.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:5cefd019c1f034e0f812f4b933a9f8973ec2df66e093e3c8370d53231479bac8 -size 6381992 diff --git a/data/ALTUM1SET/000/IMG_0245_5.tif b/data/ALTUM1SET/000/IMG_0245_5.tif deleted file mode 100644 index a60d81a4..00000000 --- a/data/ALTUM1SET/000/IMG_0245_5.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:55429bfe492a88a2e9aec9275b422e3f27cf9d9c34d6562c71c57cdca3973099 -size 6382006 diff --git a/data/ALTUM1SET/000/IMG_0245_6.tif b/data/ALTUM1SET/000/IMG_0245_6.tif deleted file mode 100644 index 735be2fb..00000000 --- a/data/ALTUM1SET/000/IMG_0245_6.tif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:a85d954314f9fa8dd6dc9175220f030d8f6ffc1c3d4fcf322b79a7fecd6b85f8 -size 44728 diff --git a/index.ipynb b/index.ipynb deleted file mode 100644 index c1ed80ce..00000000 --- a/index.ipynb +++ /dev/null @@ -1,74 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# MicaSense Image Processing Tutorials \n", - "\n", - "This repository includes tutorials and examples for processing MicaSense RedEdge and MicaSense Altum images into usable information using the python programming language. The intended audience is researchers and engineers with some software development experience that want to do their own image processing. While a number of commercial tools fully support processing RedEdge data into reflectance maps, there are a number of reasons to process your own data, including controlling the entire radiometric workflow (for academic or publication reasons), pre-processing images to be used in a non-radiometric photogrammetry suite, or processing single sets of 5 images without building a larger map.\n", - "\n", - "To get started, browse the examples below, and then head over to the [git repository](https://github.com/micasense/imageprocessing) to check out the code.\n", - "\n", - "## Tutorial Articles\n", - "\n", - "### Setup\n", - "\n", - "The setup tutorial will get you setup with the right tools to run these notebooks yourself.\n", - "\n", - "* [MicaSense Image Processing Setup](MicaSense%20Image%20Processing%20Setup.html)\n", - "\n", - "### Fundamentals\n", - "\n", - "* [Basic Radiometric Corrections](MicaSense%20Image%20Processing%20Tutorial%201.html)\n", - "* [Advanced DLS1 Processing](MicaSense%20Image%20Processing%20Tutorial%203.html)\n", - "\n", - "### micasense Library Usage\n", - "\n", - "* [MicaSense Library Introduction](MicaSense%20Image%20Processing%20Tutorial%202.html)\n", - "* [Image Class Extended Examples](Images.html)\n", - "* [Capture Class Examples](Captures.html)\n", - "* [ImageSet Examples](ImageSets.html)\n", - "* [Panels Example](Panels.html)\n", - "\n", - "### Advanced Analysis\n", - "* [Capture Alignment and Analysis Examples](Alignment.html)\n", - "* [Capture Alignment Using Rig Relatives](Alignment-RigRelatives.html)\n", - "* [Batch processing into stacked TIFFs](Batch%20Processing.html)\n", - "\n", - "### 10-Band (Dual Camera) Examples\n", - "* [10-band Capture Alignment](Alignment-10Band.html)\n", - "* [Batch Processing-10Band](Batch%20Processing-10Band.html)\n" - ] - }, - { - "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" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/micasense/image.py b/micasense/image.py index ffb94799..33c4308d 100644 --- a/micasense/image.py +++ b/micasense/image.py @@ -36,32 +36,35 @@ import micasense.metadata as metadata import micasense.dls as dls -#helper function to convert euler angles to a rotation matrix + +# helper function to convert euler angles to a rotation matrix def rotations_degrees_to_rotation_matrix(rotation_degrees): - cx = np.cos(np.deg2rad(rotation_degrees[0])) - cy = np.cos(np.deg2rad(rotation_degrees[1])) - cz = np.cos(np.deg2rad(rotation_degrees[2])) - sx = np.sin(np.deg2rad(rotation_degrees[0])) - sy = np.sin(np.deg2rad(rotation_degrees[1])) - sz = np.sin(np.deg2rad(rotation_degrees[2])) - - Rx = np.mat([ 1, 0, 0, - 0, cx,-sx, - 0, sx, cx]).reshape(3,3) - Ry = np.mat([ cy, 0, sy, - 0, 1, 0, - -sy, 0, cy]).reshape(3,3) - Rz = np.mat([ cz,-sz, 0, - sz, cz, 0, - 0, 0, 1]).reshape(3,3) - R = Rx*Ry*Rz - return R + cx = np.cos(np.deg2rad(rotation_degrees[0])) + cy = np.cos(np.deg2rad(rotation_degrees[1])) + cz = np.cos(np.deg2rad(rotation_degrees[2])) + sx = np.sin(np.deg2rad(rotation_degrees[0])) + sy = np.sin(np.deg2rad(rotation_degrees[1])) + sz = np.sin(np.deg2rad(rotation_degrees[2])) + + Rx = np.mat([1, 0, 0, + 0, cx, -sx, + 0, sx, cx]).reshape(3, 3) + Ry = np.mat([cy, 0, sy, + 0, 1, 0, + -sy, 0, cy]).reshape(3, 3) + Rz = np.mat([cz, -sz, 0, + sz, cz, 0, + 0, 0, 1]).reshape(3, 3) + R = Rx * Ry * Rz + return R + class Image(object): """ An Image is a single file taken by a RedEdge camera representing one band of multispectral information """ + def __init__(self, image_path, exiftool_obj=None): if not os.path.isfile(image_path): raise IOError("Provided path is not a file: {}".format(image_path)) @@ -72,7 +75,7 @@ def __init__(self, image_path, exiftool_obj=None): raise ValueError("Provided file path does not have a band name: {}".format(image_path)) if self.meta.band_name().upper() != 'LWIR' and not self.meta.supports_radiometric_calibration(): raise ValueError('Library requires images taken with RedEdge-(3/M/MX) camera firmware v2.1.0 or later. ' + - 'Upgrade your camera firmware to at least version 2.1.0 to use this library with RedEdge-(3/M/MX) cameras.') + 'Upgrade your camera firmware to at least version 2.1.0 to use this library with RedEdge-(3/M/MX) cameras.') self.utc_time = self.meta.utc_time() self.latitude, self.longitude, self.altitude = self.meta.position() @@ -90,8 +93,6 @@ def __init__(self, image_path, exiftool_obj=None): self.gain = self.meta.gain() self.bits_per_pixel = self.meta.bits_per_pixel() - self.vignette_center = self.meta.vignette_center() - self.vignette_polynomial = self.meta.vignette_polynomial() self.distortion_parameters = self.meta.distortion_parameters() self.principal_point = self.meta.principal_point() self.focal_plane_resolution_px_per_mm = self.meta.focal_plane_resolution_px_per_mm() @@ -101,6 +102,8 @@ def __init__(self, image_path, exiftool_obj=None): self.bandwidth = self.meta.bandwidth() self.rig_relatives = self.meta.rig_relatives() self.spectral_irradiance = self.meta.spectral_irradiance() + self.vignette_polynomial_2d = self.meta.vignette_polynomial_2d() + self.vignette_polynomial_2d_name = self.meta.vignette_polynomial_2d_name() self.auto_calibration_image = self.meta.auto_calibration_image() self.panel_albedo = self.meta.panel_albedo() @@ -108,15 +111,15 @@ def __init__(self, image_path, exiftool_obj=None): self.panel_serial = self.meta.panel_serial() if self.dls_present: - self.dls_orientation_vector = np.array([0,0,-1]) + self.dls_orientation_vector = np.array([0, 0, -1]) self.sun_vector_ned, \ self.sensor_vector_ned, \ self.sun_sensor_angle, \ self.solar_elevation, \ - self.solar_azimuth=dls.compute_sun_angle(self.location, - self.meta.dls_pose(), - self.utc_time, - self.dls_orientation_vector) + self.solar_azimuth = dls.compute_sun_angle(self.location, + self.meta.dls_pose(), + self.utc_time, + self.dls_orientation_vector) self.angular_correction = dls.fresnel(self.sun_sensor_angle) # when we have good horizontal irradiance the camera provides the solar az and el also @@ -129,37 +132,37 @@ def __init__(self, image_path, exiftool_obj=None): self.estimated_direct_vector = self.meta.estimated_direct_vector() if self.meta.horizontal_irradiance_valid(): self.horizontal_irradiance = self.meta.horizontal_irradiance() - else: + else: self.horizontal_irradiance = self.compute_horizontal_irradiance_dls2() else: - self.direct_to_diffuse_ratio = 6.0 # assumption + self.direct_to_diffuse_ratio = 6.0 # assumption self.horizontal_irradiance = self.compute_horizontal_irradiance_dls1() self.spectral_irradiance = self.meta.spectral_irradiance() - else: # no dls present or LWIR band: compute what we can, set the rest to 0 - self.dls_orientation_vector = np.array([0,0,-1]) + else: # no dls present or LWIR band: compute what we can, set the rest to 0 + self.dls_orientation_vector = np.array([0, 0, -1]) self.sun_vector_ned, \ self.sensor_vector_ned, \ self.sun_sensor_angle, \ self.solar_elevation, \ - self.solar_azimuth=dls.compute_sun_angle(self.location, - (0,0,0), - self.utc_time, - self.dls_orientation_vector) + self.solar_azimuth = dls.compute_sun_angle(self.location, + (0, 0, 0), + self.utc_time, + self.dls_orientation_vector) self.angular_correction = dls.fresnel(self.sun_sensor_angle) self.horizontal_irradiance = 0 self.scattered_irradiance = 0 self.direct_irradiance = 0 self.direct_to_diffuse_ratio = 0 - + # Internal image containers; these can use a lot of memory, clear with Image.clear_images - self.__raw_image = None # pure raw pixels - self.__intensity_image = None # black level and gain-exposure/radiometric compensated - self.__radiance_image = None # calibrated to radiance - self.__reflectance_image = None # calibrated to reflectance (0-1) + self.__raw_image = None # pure raw pixels + self.__intensity_image = None # black level and gain-exposure/radiometric compensated + self.__radiance_image = None # calibrated to radiance + self.__reflectance_image = None # calibrated to reflectance (0-1) self.__reflectance_irradiance = None - self.__undistorted_source = None # can be any of raw, intensity, radiance - self.__undistorted_image = None # current undistorted image, depdining on source + self.__undistorted_source = None # can be any of raw, intensity, radiance + self.__undistorted_image = None # current undistorted image, depdining on source # solar elevation is defined as the angle betwee the horizon and the sun, so it is 0 when the # sun is at the horizon and pi/2 when the sun is directly overhead @@ -167,25 +170,28 @@ def horizontal_irradiance_from_direct_scattered(self): return self.direct_irradiance * np.sin(self.solar_elevation) + self.scattered_irradiance def compute_horizontal_irradiance_dls1(self): - percent_diffuse = 1.0/self.direct_to_diffuse_ratio - #percent_diffuse = 5e4/(img.center_wavelength**2) + percent_diffuse = 1.0 / self.direct_to_diffuse_ratio + # percent_diffuse = 5e4/(img.center_wavelength**2) sensor_irradiance = self.spectral_irradiance / self.angular_correction # find direct irradiance in the plane normal to the sun untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(self.sun_sensor_angle)) self.direct_irradiance = untilted_direct_irr - self.scattered_irradiance = untilted_direct_irr*percent_diffuse + self.scattered_irradiance = untilted_direct_irr * percent_diffuse # compute irradiance on the ground using the solar altitude angle return self.horizontal_irradiance_from_direct_scattered() - + + def get_ready_undistorted(self): + return self.__undistorted_image + def compute_horizontal_irradiance_dls2(self): ''' Compute the proper solar elevation, solar azimuth, and horizontal irradiance for cases where the camera system did not do it correctly ''' - _,_,_, \ + _, _, _, \ self.solar_elevation, \ - self.solar_azimuth=dls.compute_sun_angle(self.location, - (0,0,0), - self.utc_time, - np.array([0,0,-1])) + self.solar_azimuth = dls.compute_sun_angle(self.location, + (0, 0, 0), + self.utc_time, + np.array([0, 0, -1])) return self.horizontal_irradiance_from_direct_scattered() def __lt__(self, other): @@ -209,31 +215,35 @@ def raw(self): import rawpy self.__raw_image = rawpy.imread(self.path).raw_image except ImportError: - self.__raw_image = cv2.imread(self.path,-1) + self.__raw_image = cv2.imread(self.path, -1) except IOError: print("Could not open image at path {}".format(self.path)) raise return self.__raw_image - - def set_raw(self,img): + + def set_raw(self, img): ''' set raw image from input img''' self.__raw_image = img.astype(np.uint16) - - def set_undistorted(self,img): + + def set_undistorted(self, img): ''' set undistorted image from input img''' self.__undistorted_image = img.astype(np.uint16) - - def set_external_rig_relatives(self,external_rig_relatives): + + def set_undistorted_custom(self, img): + ''' set undistorted image from input img''' + self.__undistorted_image = img + + def set_external_rig_relatives(self, external_rig_relatives): self.rig_translations = external_rig_relatives['rig_translations'] - #external rig relatives are in rad + # external rig relatives are in rad self.rig_relatives = [np.rad2deg(a) for a in external_rig_relatives['rig_relatives']] - px,py = external_rig_relatives['cx'],external_rig_relatives['cy'] - fx,fy = external_rig_relatives['fx'],external_rig_relatives['fy'] + px, py = external_rig_relatives['cx'], external_rig_relatives['cy'] + fx, fy = external_rig_relatives['fx'], external_rig_relatives['fy'] rx = self.focal_plane_resolution_px_per_mm[0] ry = self.focal_plane_resolution_px_per_mm[1] - self.principal_point = [px/rx,py/ry] - self.focal_length = (fx+fy)*.5/rx - #to do - set the distortion etc. + self.principal_point = [px / rx, py / ry] + self.focal_length = (fx + fy) * .5 / rx + # to do - set the distortion etc. def clear_image_data(self): ''' clear all computed images to reduce memory overhead ''' @@ -252,8 +262,8 @@ def size(self): def reflectance(self, irradiance=None, force_recompute=False): ''' Lazy-compute and return a reflectance image provided an irradiance reference ''' if self.__reflectance_image is not None \ - and force_recompute == False \ - and (self.__reflectance_irradiance == irradiance or irradiance == None): + and force_recompute == False \ + and (self.__reflectance_irradiance == irradiance or irradiance == None): return self.__reflectance_image if irradiance is None and self.band_name != 'LWIR': if self.horizontal_irradiance != 0.0: @@ -281,12 +291,12 @@ def intensity(self, force_recompute=False): _, a2, a3 = self.radiometric_cal[0], self.radiometric_cal[1], self.radiometric_cal[2] # apply image correction methods to raw image - V, x, y = self.vignette() + V, x, y = self.vignette() # self.vignette() R = 1.0 / (1.0 + a2 * y / self.exposure_time - a3 * y) L = V * R * (image_raw - self.black_level) L[L < 0] = 0 - max_raw_dn = float(2**self.bits_per_pixel) - intensity_image = L.astype(float)/(self.gain * self.exposure_time * max_raw_dn) + max_raw_dn = float(2 ** self.bits_per_pixel) + intensity_image = L.astype(float) / (self.gain * self.exposure_time * max_raw_dn) self.__intensity_image = intensity_image.T return self.__intensity_image @@ -300,7 +310,7 @@ def radiance(self, force_recompute=False): # get image dimensions image_raw = np.copy(self.raw()).T - if(self.band_name != 'LWIR'): + if (self.band_name != 'LWIR'): # get radiometric calibration factors a1, a2, a3 = self.radiometric_cal[0], self.radiometric_cal[1], self.radiometric_cal[2] # apply image correction methods to raw image @@ -308,46 +318,87 @@ def radiance(self, force_recompute=False): R = 1.0 / (1.0 + a2 * y / self.exposure_time - a3 * y) L = V * R * (image_raw - self.black_level) L[L < 0] = 0 - max_raw_dn = float(2**self.bits_per_pixel) - radiance_image = L.astype(float)/(self.gain * self.exposure_time)*a1/max_raw_dn + max_raw_dn = float(2 ** self.bits_per_pixel) + radiance_image = L.astype(float) / (self.gain * self.exposure_time) * a1 / max_raw_dn else: - L = image_raw - (273.15*100.0) # convert to C from K + L = image_raw - (273.15 * 100.0) # convert to C from K radiance_image = L.astype(float) * 0.01 self.__radiance_image = radiance_image.T return self.__radiance_image - def vignette(self): + def vignetting(self, powers_coefficients, x, y): + value = 0.0 + + for entry in powers_coefficients: + value = value + entry[2] * math.pow(x, entry[0]) * math.pow(y, entry[1]) + + return value + + def new_vignette(self): + powers_coefficients = [] + power_items = self.vignette_polynomial_2d_name + coefficient_items = self.vignette_polynomial_2d + + for i in range(0, len(power_items), 2): + powers_coefficients.append( + (int(power_items[i]), int(power_items[i + 1]), float(coefficient_items[int(i / 2)]))) + + (rows, cols) = self.raw().shape + vignette_factor = np.ones((rows, cols), dtype=np.float32) + + for y in range(0, rows): + for x in range(0, cols): + vignette_factor[y, x] = self.vignetting(powers_coefficients, x / cols, y / rows) + + x, y = np.meshgrid(np.arange(cols), np.arange(rows)) + + x = x.T + y = y.T + + # image_data = self.raw() + # image_data = image_data / vignette_factor + final_vignette = (1 / vignette_factor).T + + return final_vignette, x, y + + def vignette(self): # 2D polynomial ''' Get a numpy array which defines the value to multiply each pixel by to correct for optical vignetting effects. Note: this array is transposed from normal image orientation and comes as part of a three-tuple, the other parts of which are also used by the radiance method. ''' - # get vignette center - vignette_center_x, vignette_center_y = self.vignette_center - - # get a copy of the vignette polynomial because we want to modify it here - v_poly_list = list(self.vignette_polynomial) - # reverse list and append 1., so that we can call with numpy polyval - v_poly_list.reverse() - v_poly_list.append(1.) - v_polynomial = np.array(v_poly_list) + # get vignette 2d polynomial terms + dim_poly = np.max(self.vignette_polynomial_2d_name) + 1 + v_poly_mat = np.zeros(shape=(dim_poly, dim_poly)) + for i in range(0, len(self.vignette_polynomial_2d)): + x_m = self.vignette_polynomial_2d_name[2 * i] + y_m = self.vignette_polynomial_2d_name[2 * i + 1] + v_poly_mat[x_m, y_m] = self.vignette_polynomial_2d[i] # perform vignette correction # get coordinate grid across image, seem swapped because of transposed vignette + x_dim, y_dim = self.raw().shape[1], self.raw().shape[0] x, y = np.meshgrid(np.arange(x_dim), np.arange(y_dim)) - #meshgrid returns transposed arrays + # meshgrid returns transposed arrays x = x.T y = y.T - # compute matrix of distances from image center - r = np.hypot((x-vignette_center_x), (y-vignette_center_y)) + w_norm = x / x_dim + h_norm = y / y_dim + + vignette_mat = np.zeros(shape=(x_dim, y_dim)) + + for i in range(v_poly_mat.shape[0]): + for j in range(v_poly_mat.shape[1]): + vignette_mat += v_poly_mat[i, j] * (w_norm ** i) * (h_norm ** j) # compute the vignette polynomial for each distance - we divide by the polynomial so that the # corrected image is image_corrected = image_original * vignetteCorrection - vignette = 1./np.polyval(v_polynomial, r) + vignette = 1. / vignette_mat + return vignette, x, y def undistorted_radiance(self, force_recompute=False): @@ -360,7 +411,7 @@ def plottable_vignette(self): return self.vignette()[0].T def cv2_distortion_coeff(self): - #dist_coeffs = np.array(k[0],k[1],p[0],p[1],k[2]]) + # dist_coeffs = np.array(k[0],k[1],p[0],p[1],k[2]]) return np.array(self.distortion_parameters)[[0, 1, 3, 4, 2]] # values in pp are in [mm], rescale to pixels @@ -384,10 +435,10 @@ def cv2_camera_matrix(self): return cam_mat def rig_xy_offset_in_px(self): - pixel_pitch_mm_x = 1.0/self.focal_plane_resolution_px_per_mm[0] - pixel_pitch_mm_y = 1.0/self.focal_plane_resolution_px_per_mm[1] - px_fov_x = 2.0 * math.atan2(pixel_pitch_mm_x/2.0, self.focal_length) - px_fov_y = 2.0 * math.atan2(pixel_pitch_mm_y/2.0, self.focal_length) + pixel_pitch_mm_x = 1.0 / self.focal_plane_resolution_px_per_mm[0] + pixel_pitch_mm_y = 1.0 / self.focal_plane_resolution_px_per_mm[1] + px_fov_x = 2.0 * math.atan2(pixel_pitch_mm_x / 2.0, self.focal_length) + px_fov_y = 2.0 * math.atan2(pixel_pitch_mm_y / 2.0, self.focal_length) t_x = math.radians(self.rig_relatives[0]) / px_fov_x t_y = math.radians(self.rig_relatives[1]) / px_fov_y return (t_x, t_y) @@ -406,11 +457,11 @@ def undistorted(self, image): self.size(), 1) map1, map2 = cv2.initUndistortRectifyMap(self.cv2_camera_matrix(), - self.cv2_distortion_coeff(), - np.eye(3), - new_cam_mat, - self.size(), - cv2.CV_32F) # cv2.CV_32F for 32 bit floats + self.cv2_distortion_coeff(), + np.eye(3), + new_cam_mat, + self.size(), + cv2.CV_32F) # cv2.CV_32F for 32 bit floats # compute the undistorted 16 bit image self.__undistorted_image = cv2.remap(image, map1, map2, cv2.INTER_LINEAR) return self.__undistorted_image @@ -427,7 +478,6 @@ def plot_intensity(self, title=None, figsize=None): title = '{} Band {} Intensity (DN*sec)'.format(self.band_name, self.band_index) return plotutils.plotwithcolorbar(self.intensity(), title=title, figsize=figsize) - def plot_radiance(self, title=None, figsize=None): ''' Create a single plot of the image converted to radiance ''' if title is None: @@ -446,41 +496,42 @@ def plot_undistorted_radiance(self, title=None, figsize=None): title = '{} Band {} Undistorted Radiance'.format(self.band_name, self.band_index) return plotutils.plotwithcolorbar(self.undistorted(self.radiance()), title=title, figsize=figsize) - def plot_all(self, figsize=(13,10)): + def plot_all(self, figsize=(13, 10)): plots = [self.raw(), self.plottable_vignette(), self.radiance(), self.undistorted(self.radiance())] plot_types = ['Raw', 'Vignette', 'Radiance', 'Undistorted Radiance'] titles = ['{} Band {} {}'.format(str(self.band_name), str(self.band_index), tpe) - for tpe in plot_types] + for tpe in plot_types] plotutils.subplotwithcolorbar(2, 2, plots, titles, figsize=figsize) - #get the homography that maps from this image to the reference image - def get_homography(self,ref,R=None,T=None): - #if we have externally supplied rotations/translations for the rig use these + # get the homography that maps from this image to the reference image + + def get_homography(self, ref, R=None, T=None): + # if we have externally supplied rotations/translations for the rig use these # otherwise use the rig-relatives intrinsic to the image if R is None: R = rotations_degrees_to_rotation_matrix(self.rig_relatives) if T is None: - T =np.zeros(3) + T = np.zeros(3) R_ref = rotations_degrees_to_rotation_matrix(ref.rig_relatives) - A = np.zeros((4,4)) - A[0:3,0:3]=np.dot(R_ref.T,R) - A[0:3,3]=T - A[3,3]=1. + A = np.zeros((4, 4)) + A[0:3, 0:3] = np.dot(R_ref.T, R) + A[0:3, 3] = T + A[3, 3] = 1. C, _ = cv2.getOptimalNewCameraMatrix(self.cv2_camera_matrix(), self.cv2_distortion_coeff(), - self.size(),1) + self.size(), 1) Cr, _ = cv2.getOptimalNewCameraMatrix(ref.cv2_camera_matrix(), ref.cv2_distortion_coeff(), - ref.size(),1) - CC = np.zeros((4,4)) - CC[0:3,0:3] = C - CC[3,3]=1. - CCr = np.zeros((4,4)) - CCr[0:3,0:3] = Cr - CCr[3,3]=1. - - B = np.array(np.dot(CCr,np.dot(A,np.linalg.inv(CC)))) - B[:,2]=B[:,2]-B[:,3] - B = B[0:3,0:3] - B = B/B[2,2] + ref.size(), 1) + CC = np.zeros((4, 4)) + CC[0:3, 0:3] = C + CC[3, 3] = 1. + CCr = np.zeros((4, 4)) + CCr[0:3, 0:3] = Cr + CCr[3, 3] = 1. + + B = np.array(np.dot(CCr, np.dot(A, np.linalg.inv(CC)))) + B[:, 2] = B[:, 2] - B[:, 3] + B = B[0:3, 0:3] + B = B / B[2, 2] return np.array(B) diff --git a/micasense/imageutils.py b/micasense/imageutils.py index 643c365b..04ea0fd8 100644 --- a/micasense/imageutils.py +++ b/micasense/imageutils.py @@ -185,7 +185,8 @@ def default_warp_matrix(warp_mode): else: return np.array([[1,0,0],[0,1,0]], dtype=np.float32) -def align_capture(capture, ref_index=1, warp_mode=cv2.MOTION_HOMOGRAPHY, max_iterations=2500, epsilon_threshold=1e-9, multithreaded=True, debug=False, pyramid_levels = None): + +def align_capture_original(capture, ref_index=1, warp_mode=cv2.MOTION_HOMOGRAPHY, max_iterations=2500, epsilon_threshold=1e-9, multithreaded=True, debug=False, pyramid_levels = None): '''Align images in a capture using openCV MOTION_TRANSLATION sets a translational motion model; warpMatrix is 2x3 with the first 2x2 part being the unity matrix and the rest two parameters being estimated. MOTION_EUCLIDEAN sets a Euclidean (rigid) transformation as motion model; three parameters are estimated; warpMatrix is 2x3. @@ -256,6 +257,238 @@ def align_capture(capture, ref_index=1, warp_mode=cv2.MOTION_HOMOGRAPHY, max_ite warp_matrices.append(capture.get_warp_matrices(ref_index)[-1]) return warp_matrices, alignment_pairs + + + +def align_capture(capture, ref_index=1, warp_modes=[], max_iterations=2500, epsilon_threshold=1e-9, multithreaded=True, debug=False, pyramid_levels = None, img_type='reflectance'): + '''Align images in a capture using openCV # cv2.MOTION_HOMOGRAPHY + MOTION_TRANSLATION sets a translational motion model; warpMatrix is 2x3 with the first 2x2 part being the unity matrix and the rest two parameters being estimated. + MOTION_EUCLIDEAN sets a Euclidean (rigid) transformation as motion model; three parameters are estimated; warpMatrix is 2x3. + MOTION_AFFINE sets an affine motion model (DEFAULT); six parameters are estimated; warpMatrix is 2x3. + MOTION_HOMOGRAPHY sets a homography as a motion model; eight parameters are estimated;`warpMatrix` is 3x3. + best results will be AFFINE and HOMOGRAPHY, at the expense of speed + ''' + # Match other bands to this reference image (index into capture.images[]) + ref_img = capture.images[ref_index].undistorted(capture.images[ref_index].radiance()).astype('float32') + + if capture.has_rig_relatives(): + warp_matrices_init = capture.get_warp_matrices(ref_index=ref_index) + else: + warp_matrices_init = [default_warp_matrix(warp_modes[0])]*len(capture.images) + + alignment_pairs = [] + for i,img in enumerate(capture.images): + if img.rig_relatives is not None: + translations = img.rig_xy_offset_in_px() + else: + translations = (0,0) + if img.band_name != 'LWIR': + alignment_pairs.append({'warp_mode': warp_modes[0], + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index':ref_index, + 'ref_image': ref_img, + 'match_index':i, + 'match_image':img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'warp_matrix_init': np.array(warp_matrices_init[i], dtype=np.float32), + 'debug': debug, + 'pyramid_levels': pyramid_levels}) + warp_matrices = [None]*len(alignment_pairs) + + #required to work across linux/mac/windows, see https://stackoverflow.com/questions/47852237 + if multithreaded and multiprocessing.get_start_method() != 'spawn': + try: + multiprocessing.set_start_method('spawn',force=True) + except ValueError: + multithreaded = False + + if(multithreaded): + pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + for _, mat in enumerate(pool.imap_unordered(align, alignment_pairs)): + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + pool.close() + pool.join() + else: + # Single-threaded alternative + for pair in alignment_pairs: + mat = align(pair) + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + + if capture.images[-1].band_name == 'LWIR': + img = capture.images[-1] + alignment_pairs.append({'warp_mode': warp_modes[0], + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index':ref_index, + 'ref_image': ref_img, + 'match_index':img.band_index, + 'match_image':img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'debug': debug}) + warp_matrices.append(capture.get_warp_matrices(ref_index)[-1]) + + # Create aligned stack + # cropped_dimensions, edges = find_crop_bounds(capture, warp_matrices, warp_mode=warp_modes[0]) + im_aligned = aligned_capture_without_image_cropping(capture, warp_matrices, warp_modes[0], img_type=img_type) + + # Another alignment + alignment_pairs = [] + for i, img in enumerate(capture.images): + # if img.rig_relatives is not None: + # translations = img.rig_xy_offset_in_px() + # else: + translations = (0, 0) + if img.band_name != 'LWIR': + alignment_pairs.append({'warp_mode': warp_modes[1], + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index': ref_index, + 'ref_image': ref_img, + 'match_index': i, + 'match_image': im_aligned[:, :, i], # img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'warp_matrix_init': np.array(warp_matrices_init[i], dtype=np.float32), + 'debug': debug, + 'pyramid_levels': pyramid_levels}) + + warp_matrices = [None] * len(alignment_pairs) + + # required to work across linux/mac/windows, see https://stackoverflow.com/questions/47852237 + if multithreaded and multiprocessing.get_start_method() != 'spawn': + try: + multiprocessing.set_start_method('spawn', force=True) + except ValueError: + multithreaded = False + + if (multithreaded): + pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + for _, mat in enumerate(pool.imap_unordered(align, alignment_pairs)): + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + pool.close() + pool.join() + else: + # Single-threaded alternative + for pair in alignment_pairs: + mat = align(pair) + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + + if capture.images[-1].band_name == 'LWIR': + img = capture.images[-1] + alignment_pairs.append({'warp_mode': warp_modes[1], + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index': ref_index, + 'ref_image': ref_img, + 'match_index': img.band_index, + 'match_image': img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'debug': debug}) + warp_matrices.append(capture.get_warp_matrices(ref_index)[-1]) + + + return warp_matrices, alignment_pairs + + +def align_capture_custom(capture, ref_index=1, warp_mode=cv2.MOTION_HOMOGRAPHY, max_iterations=2500, epsilon_threshold=1e-9, multithreaded=True, debug=False, pyramid_levels = None, ignore_rig = False): + '''Align images in a capture using openCV + MOTION_TRANSLATION sets a translational motion model; warpMatrix is 2x3 with the first 2x2 part being the unity matrix and the rest two parameters being estimated. + MOTION_EUCLIDEAN sets a Euclidean (rigid) transformation as motion model; three parameters are estimated; warpMatrix is 2x3. + MOTION_AFFINE sets an affine motion model (DEFAULT); six parameters are estimated; warpMatrix is 2x3. + MOTION_HOMOGRAPHY sets a homography as a motion model; eight parameters are estimated;warpMatrix is 3x3. + best results will be AFFINE and HOMOGRAPHY, at the expense of speed + ''' + # Match other bands to this reference image (index into capture.images[]) + ref_img = capture.images[ref_index].undistorted(capture.images[ref_index].radiance()).astype('float32') + if not ignore_rig and capture.has_rig_relatives(): + warp_matrices_init = capture.get_warp_matrices(ref_index=ref_index) + else: + warp_matrices_init = [default_warp_matrix(warp_mode)]*len(capture.images) + alignment_pairs = [] + for i,img in enumerate(capture.images): + if not ignore_rig and img.rig_relatives is not None: + translations = img.rig_xy_offset_in_px() + else: + translations = (0,0) + if img.band_name != 'LWIR': + alignment_pairs.append({'warp_mode': warp_mode, + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index':ref_index, + 'ref_image': ref_img, + 'match_index':i, + 'match_image':img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'warp_matrix_init': np.array(warp_matrices_init[i], dtype=np.float32), + 'debug': debug, + 'pyramid_levels': pyramid_levels}) + + warp_matrices = [None]*len(alignment_pairs) + + #required to work across linux/mac/windows, see https://stackoverflow.com/questions/47852237 + if multithreaded and multiprocessing.get_start_method() != 'spawn': + try: + multiprocessing.set_start_method('spawn',force=True) + except ValueError: + multithreaded = False + + if(multithreaded): + pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) + for _,mat in enumerate(pool.imap_unordered(align, alignment_pairs)): + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + pool.close() + pool.join() + else: + # Single-threaded alternative + for pair in alignment_pairs: + mat = align(pair) + warp_matrices[mat['match_index']] = mat['warp_matrix'] + print("Finished aligning band {}".format(mat['match_index'])) + if capture.images[-1].band_name == 'LWIR': + img = capture.images[-1] + alignment_pairs.append({'warp_mode': warp_mode, + 'max_iterations': max_iterations, + 'epsilon_threshold': epsilon_threshold, + 'ref_index':ref_index, + 'ref_image': ref_img, + 'match_index':img.band_index, + 'match_image':img.undistorted(img.radiance()).astype('float32'), + 'translations': translations, + 'debug': debug}) + warp_matrices.append(capture.get_warp_matrices(ref_index)[-1]) + + return warp_matrices, alignment_pairs + + +#apply homography to create an aligned stack +def aligned_capture_without_image_cropping(capture, warp_matrices, warp_mode, img_type='reflectance', interpolation_mode=cv2.INTER_LANCZOS4): + width, height = capture.images[0].size() + + im_aligned = np.zeros((height,width,len(warp_matrices)), dtype=np.float32 ) + + for i in range(0,len(warp_matrices)): + if img_type == 'reflectance': + img = capture.images[i].undistorted_reflectance() + else: + img = capture.images[i].undistorted_radiance() + + if warp_mode != cv2.MOTION_HOMOGRAPHY: + im_aligned[:,:,i] = cv2.warpAffine(img, + warp_matrices[i], + (width,height), + flags=interpolation_mode + cv2.WARP_INVERSE_MAP) + else: + im_aligned[:,:,i] = cv2.warpPerspective(img, + warp_matrices[i], + (width,height), + flags=interpolation_mode + cv2.WARP_INVERSE_MAP) + return im_aligned + #apply homography to create an aligned stack def aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, match_index, img_type = 'reflectance',interpolation_mode=cv2.INTER_LANCZOS4): width, height = capture.images[0].size() diff --git a/micasense/metadata.py b/micasense/metadata.py index fbffb249..edb34fd7 100644 --- a/micasense/metadata.py +++ b/micasense/metadata.py @@ -30,10 +30,14 @@ import pytz import os import math +import random class Metadata(object): ''' Container for Micasense image metadata''' def __init__(self, filename, exiftoolPath=None, exiftool_obj=None): + self.calculated_vignette_center = 0 + self.vignette_center_x = 0 + self.vignette_center_y = 0 if exiftool_obj is not None: self.exif = exiftool_obj.get_metadata(filename) return @@ -228,15 +232,33 @@ def bits_per_pixel(self): ''' get the number of bits per pixel, which defines the maximum digital number value in the image ''' return self.get_item('EXIF:BitsPerSample') - def vignette_center(self): - ''' get the vignette center in X and Y image coordinates''' - nelem = self.size('XMP:VignettingCenter') - return [float(self.get_item('XMP:VignettingCenter', i)) for i in range(nelem)] - - def vignette_polynomial(self): + # def _calculate_vignette_center(self): + # if not self.calculated_vignette_center: + # self.calculated_vignette_center = 1 + # image_width = self.get_item('EXIF:ImageWidth') + # image_height = self.get_item('EXIF:ImageHeight') + # random_gaussian_w = random.gauss(2, 0.1) + # random_gaussian_h = random.gauss(2, 0.1) + # + # self.vignette_center_x = image_width / 2 # random_gaussian_w + # self.vignette_center_y = image_height / 2 # random_gaussian_h + # + # return [self.vignette_center_x, self.vignette_center_y] + # else: + # return [self.vignette_center_x, self.vignette_center_y] + # + # def vignette_center(self): + # ''' get the vignette center in X and Y image coordinates''' + # return self._calculate_vignette_center() + + def vignette_polynomial_2d_name(self): + poly_2d_name = self.get_item('XMP:VignettingPolynomial2DName') + return [ int(poly_v) for poly_v in poly_2d_name.split(',')] + + def vignette_polynomial_2d(self): ''' get the radial vignette polynomial in the order it's defined in the metadata''' - nelem = self.size('XMP:VignettingPolynomial') - return [float(self.get_item('XMP:VignettingPolynomial', i)) for i in range(nelem)] + poly_2d = self.get_item('XMP:VignettingPolynomial2D') + return [ float(poly_v) for poly_v in poly_2d.split(',')] def distortion_parameters(self): nelem = self.size('XMP:PerspectiveDistortion') @@ -295,7 +317,7 @@ def horizontal_irradiance_valid(self): version_string = self.firmware_version().strip('v') if self.camera_model() == "Altum": good_version = "1.2.3" - elif self.camera_model() == 'RedEdge' or self.camera_model() == 'RedEdge-M': + elif self.camera_model() == 'RedEdge' or self.camera_model() == 'RedEdge-M' or self.camera_model() == 'RedEdge-P': good_version = "5.1.7" else: raise ValueError("Camera model is required to be RedEdge or Altum, not {} ".format(self.camera_model())) diff --git a/micasense/panel.py b/micasense/panel.py index 2c314afe..4b8a63d2 100644 --- a/micasense/panel.py +++ b/micasense/panel.py @@ -194,7 +194,7 @@ def panel_corners(self): reference_qr_pts = np.asarray([[-p, p], [p, p], [p, -p], [-p, -p]], dtype=np.float32) bounds = [] costs = [] - for rotation in range(0,4): + for rotation in [3]: qr_points = np.roll(reference_qr_pts, rotation, axis=0) src = np.asarray([tuple(row) for row in qr_points[:]], np.float32) diff --git a/micasense/utils.py b/micasense/utils.py index 93006b8b..c6c203be 100644 --- a/micasense/utils.py +++ b/micasense/utils.py @@ -78,12 +78,12 @@ def raw_image_to_radiance(meta, imageRaw): def vignette_map(meta, xDim, yDim): # get vignette center - xVignette = float(meta.get_item('XMP:VignettingCenter', 0)) - yVignette = float(meta.get_item('XMP:VignettingCenter', 1)) + xVignette = meta.vignette_center()[0] # float(meta.get_item('XMP:VignettingCenter', 0)) + yVignette = meta.vignette_center()[1] # float(meta.get_item('XMP:VignettingCenter', 1)) # get vignette polynomial - NvignettePoly = meta.size('XMP:VignettingPolynomial') - vignettePolyList = [float(meta.get_item('XMP:VignettingPolynomial', i)) for i in range(NvignettePoly)] + NvignettePoly = meta.size('XMP:VignettingPolynomial2D') + vignettePolyList = [float(meta.get_item('XMP:VignettingPolynomial2D', i)) for i in range(NvignettePoly)] # reverse list and append 1., so that we can call with numpy polyval vignettePolyList.reverse()