From cc86cc31775f36e2de97ff218910ba0e30122a9a Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 10:53:12 +0000 Subject: [PATCH 1/9] Fix SIFT alignment and use spawn pools for rawpy safety. Correct reference-band handling and calibrated warp fallback in SIFT_align_capture, add spawn_pool() for library multiprocessing, declare rawpy, and add pipeline compatibility plan plus tests. --- CHANGELOG.md | 8 ++++++++ micasense/capture.py | 18 +++++++++++------- micasense/imageset.py | 15 ++++++++------- micasense/imageutils.py | 17 +++++------------ micasense/mp_config.py | 10 ++++++++++ pyproject.toml | 1 + tests/test_sift_align.py | 33 +++++++++++++++++++++++++++++++++ 7 files changed, 76 insertions(+), 26 deletions(-) create mode 100644 micasense/mp_config.py create mode 100644 tests/test_sift_align.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7156a90..57db22c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * repr functions to classes so that debugging is easier * delete old config files * fixed jupyter notebooks +* `micasense.mp_config.spawn_pool()` using a per-pool spawn context (rawpy/OpenMP safe on Linux) + +### Fixed + +* `SIFT_align_capture`: assign reference SIFT image when reference and target bands share shape (fixes `NameError` when `ref` is a multispectral band) +* `SIFT_align_capture`: fall back to calibrated warp matrices when SIFT match count is low (e.g. LWIR) instead of raising +* `ImageSet.save_stacks` and `imageutils` alignment pools use spawn context (not global `set_start_method`) +* declare `rawpy` as a runtime dependency ## [0.1.1] - 2025-12-28 diff --git a/micasense/capture.py b/micasense/capture.py index a4cda72..1326aa9 100644 --- a/micasense/capture.py +++ b/micasense/capture.py @@ -1167,12 +1167,10 @@ def SIFT_align_capture( # use the calibrated warp matrices to verify keypoints warp_matrices_calibrated = self.get_warp_matrices(ref_index=ref) - ref_image_SIFT = self.images[ref].undistorted(self.images[ref].raw()) - if rest_shape != ref_shape: - ref_image_SIFT = resize(ref_image_SIFT, rest_shape) - ref_image_SIFT = (ref_image_SIFT / ref_image_SIFT.max() * 65535).astype( - np.uint16 - ) + ref_img = self.images[ref].undistorted(self.images[ref].raw()) + if ref_shape != rest_shape: + ref_img = resize(ref_img, rest_shape) + ref_image_SIFT = (ref_img / ref_img.max() * 65535).astype(np.uint16) descriptor_extractor.detect_and_extract(ref_image_SIFT) keypoints_ref = descriptor_extractor.keypoints @@ -1283,7 +1281,13 @@ def SIFT_align_capture( # most of the time this will occur for the thermal image, as we have a hard time # finding a good matches between panchro & thermal in most cases else: - raise Exception(f"no match for index {ix}") + P = ProjectiveTransform(matrix=warp_matrices_calibrated[ix]) + logger.warning( + "Insufficient SIFT matches for band index %s (%s); using calibrated warp matrix.", + ix, + self.images[ix].band_name, + ) + kpi, kpr = filtered_kpi, filtered_kpr models.append(P) kp_image.append(kpi) kp_ref.append(kpr) diff --git a/micasense/imageset.py b/micasense/imageset.py index a8e5ecc..36ce038 100644 --- a/micasense/imageset.py +++ b/micasense/imageset.py @@ -29,9 +29,10 @@ """ import fnmatch -import multiprocessing import os +from micasense.mp_config import spawn_pool + import exiftool import micasense.capture as capture @@ -170,12 +171,12 @@ def save_stacks( ) if multiprocess: - pool = multiprocessing.Pool(processes=multiprocessing.cpu_count()) - for i, _ in enumerate(pool.imap_unordered(save_capture, save_params_list)): - if progress_callback is not None: - progress_callback(float(i) / float(len(save_params_list))) - pool.close() - pool.join() + with spawn_pool() as pool: + for i, _ in enumerate( + pool.imap_unordered(save_capture, save_params_list) + ): + if progress_callback is not None: + progress_callback(float(i) / float(len(save_params_list))) else: for params in save_params_list: save_capture(params) diff --git a/micasense/imageutils.py b/micasense/imageutils.py index 17dbd29..f1607a4 100644 --- a/micasense/imageutils.py +++ b/micasense/imageutils.py @@ -387,20 +387,13 @@ def align_capture( ) 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 + from micasense.mp_config import spawn_pool 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"] - logger.info("Finished aligning band %s", mat["match_index"]) - pool.close() - pool.join() + with spawn_pool() as pool: + for _, mat in enumerate(pool.imap_unordered(align, alignment_pairs)): + warp_matrices[mat["match_index"]] = mat["warp_matrix"] + logger.info("Finished aligning band %s", mat["match_index"]) else: # Single-threaded alternative for pair in alignment_pairs: diff --git a/micasense/mp_config.py b/micasense/mp_config.py new file mode 100644 index 0000000..16ca8ae --- /dev/null +++ b/micasense/mp_config.py @@ -0,0 +1,10 @@ +"""Multiprocessing helpers for rawpy/OpenMP safety.""" + +import multiprocessing as mp +from multiprocessing.pool import Pool +from typing import Optional + + +def spawn_pool(processes: Optional[int] = None) -> Pool: + """Process pool using spawn (safe with rawpy/libraw OpenMP on Linux).""" + return mp.get_context("spawn").Pool(processes=processes or mp.cpu_count()) diff --git a/pyproject.toml b/pyproject.toml index 058f969..b506f16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,7 @@ dependencies = [ "requests", "numpy", "opencv-python", + "rawpy", "gdal", "pysolar", "matplotlib", diff --git a/tests/test_sift_align.py b/tests/test_sift_align.py new file mode 100644 index 0000000..c2e4b4b --- /dev/null +++ b/tests/test_sift_align.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python +"""Tests for SIFT_align_capture and multiprocessing spawn configuration.""" + +import numpy as np +import pytest + +from micasense.capture import Capture +from micasense.mp_config import spawn_pool + + +def test_sift_align_smoke(non_panel_rededge_capture): + matrices = non_panel_rededge_capture.SIFT_align_capture(ref=0, min_matches=10) + assert len(matrices) == len(non_panel_rededge_capture.images) + assert np.allclose(matrices[0], np.eye(3)) + + +def test_sift_align_same_shape_ref(non_panel_rededge_capture): + """MS band as ref when all bands share resolution must not raise NameError.""" + matrices = non_panel_rededge_capture.SIFT_align_capture(ref=2, min_matches=10) + assert len(matrices) == len(non_panel_rededge_capture.images) + assert np.allclose(matrices[2], np.eye(3)) + + +def test_sift_align_lwir_fallback(non_panel_altum_capture): + """Altum LWIR should fall back to calibrated warp instead of aborting.""" + matrices = non_panel_altum_capture.SIFT_align_capture(ref=0, min_matches=10) + assert len(matrices) == len(non_panel_altum_capture.images) + assert np.allclose(matrices[0], np.eye(3)) + + +def test_spawn_pool_uses_spawn_context(): + with spawn_pool(1) as pool: + assert pool._ctx.get_start_method() == "spawn" From c3ecce45df37ec1e7b4fe70c720451edd01c8675 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 15:04:59 +0000 Subject: [PATCH 2/9] Add warp matrix save/load API and remove dead SIFT cache field. Introduce micasense.warp_io for standard .npy warp matrix I/O and drop unused Capture.__sift_warp_matrices. --- CHANGELOG.md | 2 ++ micasense/capture.py | 1 - micasense/warp_io.py | 52 +++++++++++++++++++++++++++++++++++++++++++ tests/test_warp_io.py | 43 +++++++++++++++++++++++++++++++++++ 4 files changed, 97 insertions(+), 1 deletion(-) create mode 100644 micasense/warp_io.py create mode 100644 tests/test_warp_io.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 57db22c..6c27f7f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * delete old config files * fixed jupyter notebooks * `micasense.mp_config.spawn_pool()` using a per-pool spawn context (rawpy/OpenMP safe on Linux) +* `micasense.warp_io` — `save_warp_matrices` / `load_warp_matrices` for `.npy` warp matrix I/O ### Fixed @@ -20,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * `SIFT_align_capture`: fall back to calibrated warp matrices when SIFT match count is low (e.g. LWIR) instead of raising * `ImageSet.save_stacks` and `imageutils` alignment pools use spawn context (not global `set_start_method`) * declare `rawpy` as a runtime dependency +* remove unused `Capture.__sift_warp_matrices` ## [0.1.1] - 2025-12-28 diff --git a/micasense/capture.py b/micasense/capture.py index 1326aa9..0f5bdb6 100644 --- a/micasense/capture.py +++ b/micasense/capture.py @@ -104,7 +104,6 @@ def __init__(self, images, panel_corners=None): self.__aligned_capture = None self.__aligned_radiometric_pan_sharpened_capture = None - self.__sift_warp_matrices = None def set_panel_corners(self, panel_corners): """ diff --git a/micasense/warp_io.py b/micasense/warp_io.py new file mode 100644 index 0000000..c40d77c --- /dev/null +++ b/micasense/warp_io.py @@ -0,0 +1,52 @@ +"""Load and save band alignment warp matrices (.npy).""" + +from pathlib import Path +from typing import List, Union + +import numpy as np +from skimage.transform import ProjectiveTransform + +MatrixLike = Union[np.ndarray, ProjectiveTransform] + + +def _matrix_to_array(matrix: MatrixLike) -> np.ndarray: + if isinstance(matrix, ProjectiveTransform): + return np.asarray(matrix.params, dtype=np.float64) + return np.asarray(matrix, dtype=np.float64) + + +def warp_matrices_to_arrays(matrices: List[MatrixLike]) -> np.ndarray: + """Convert warp matrices to an object array of 3×3 float64 (on-disk format).""" + return np.array([_matrix_to_array(m) for m in matrices], dtype=object) + + +def arrays_to_warp_matrices( + arrays: np.ndarray, *, as_projective: bool = True +) -> list: + """ + Convert a loaded object array back to warp matrices. + + :param as_projective: If True, return ``ProjectiveTransform`` (SIFT / skimage). + If False, return ``ndarray`` (OpenCV alignment path). + """ + result = [] + for item in arrays: + arr = np.asarray(item) + if as_projective: + result.append(ProjectiveTransform(matrix=arr.astype(np.float64))) + else: + result.append(arr.astype(np.float32)) + return result + + +def save_warp_matrices(path: Union[Path, str], matrices: List[MatrixLike]) -> None: + """Save warp matrices in the standard ``.npy`` format (object array, pickle).""" + np.save(Path(path), warp_matrices_to_arrays(matrices), allow_pickle=True) + + +def load_warp_matrices( + path: Union[Path, str], *, as_projective: bool = True +) -> list: + """Load warp matrices written by :func:`save_warp_matrices`.""" + arrays = np.load(Path(path), allow_pickle=True) + return arrays_to_warp_matrices(arrays, as_projective=as_projective) diff --git a/tests/test_warp_io.py b/tests/test_warp_io.py new file mode 100644 index 0000000..7511cff --- /dev/null +++ b/tests/test_warp_io.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +"""Tests for warp matrix save/load.""" + +import numpy as np +from skimage.transform import ProjectiveTransform + +from micasense.warp_io import ( + arrays_to_warp_matrices, + load_warp_matrices, + save_warp_matrices, + warp_matrices_to_arrays, +) + + +def test_warp_matrices_round_trip_ndarray(panel_altum_capture, tmp_path): + matrices = panel_altum_capture.get_warp_matrices() + path = tmp_path / "warp.npy" + save_warp_matrices(path, matrices) + loaded = load_warp_matrices(path) + for orig, transform in zip(matrices, loaded): + np.testing.assert_allclose(orig, transform.params, atol=1e-6) + + +def test_warp_matrices_round_trip_projective(tmp_path): + matrices = [ + ProjectiveTransform(matrix=np.eye(3)), + ProjectiveTransform(matrix=np.diag([1.0, 1.0, 1.0])), + ] + path = tmp_path / "warp_sift.npy" + save_warp_matrices(path, matrices) + loaded = load_warp_matrices(path) + for orig, loaded_t in zip(matrices, loaded): + np.testing.assert_allclose(orig.params, loaded_t.params, atol=1e-6) + + +def test_load_as_arrays(tmp_path): + arrays = warp_matrices_to_arrays([np.eye(3), np.eye(3) * 2]) + path = tmp_path / "warp_cv.npy" + np.save(path, arrays, allow_pickle=True) + loaded = load_warp_matrices(path, as_projective=False) + assert len(loaded) == 2 + assert all(isinstance(m, np.ndarray) for m in loaded) + np.testing.assert_allclose(loaded[0], np.eye(3), atol=1e-6) From 0ec78ab4f75d0f2384094549b7c8c8bcbde9bbf0 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 15:07:23 +0000 Subject: [PATCH 3/9] Use warp_io for warp matrix load/save in tutorials and batch script. Replace duplicated np.load/np.save loops in Alignment v2, Batch Processing v2, and batch_processing_script with micasense.warp_io helpers. Remove hardcoded Mapbox token from Batch Processing v2. --- Alignment v2.ipynb | 1894 ++++++++++++++++++------------------ Batch Processing v2.ipynb | 759 +++++++-------- CHANGELOG.md | 1 + batch_processing_script.py | 19 +- 4 files changed, 1323 insertions(+), 1350 deletions(-) diff --git a/Alignment v2.ipynb b/Alignment v2.ipynb index 41c374c..88676ac 100644 --- a/Alignment v2.ipynb +++ b/Alignment v2.ipynb @@ -1,957 +1,941 @@ { - "cells": [ - { - "cell_type": "markdown", - "id": "573422ea", - "metadata": {}, - "source": [ - "# Active Image Alignment\n", - "For 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 RedEdge (3,M,MX) versions above 3.2.0 which include the \"RigRelatives\" tags, and all RedEdge-P/Altum/Altum-PT imagery. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.\n", - "\n", - "# Opening Images\n", - "As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize all the bands of a MicaSense 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 the `autoreload` extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time." - ] + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Active Image Alignment\n", + "For 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 RedEdge (3,M,MX) versions above 3.2.0 which include the \"RigRelatives\" tags, and all RedEdge-P/Altum/Altum-PT imagery. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.\n", + "\n", + "# Opening Images\n", + "As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize all the bands of a MicaSense 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 the `autoreload` extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time." + ], + "id": "573422ea" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ], + "execution_count": null, + "outputs": [], + "id": "b09cd3b1" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import micasense.capture as capture\n", + "\n", + "%matplotlib inline\n", + "from pathlib import Path\n", + "import matplotlib.pyplot as plt\n", + "\n", + "plt.rcParams[\"figure.facecolor\"] = \"w\"\n", + "\n", + "panelNames = None\n", + "\n", + "# set your image paths here. See more here: https://docs.python.org/3/library/pathlib.html\n", + "# if using Windows, you need to an an \"r\" to the path like this: Path(r\"C:\\Files\")\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-MX-DUAL\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0431_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-MX\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0020_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0001_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/ALTUM\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0021_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-P\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0011_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "imagePath = Path(\"./data/ALTUM-PT\")\n", + "\n", + "# these will return lists of image paths as strings\n", + "imageNames = list(imagePath.glob(\"IMG_0010_*.tif\"))\n", + "imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "panelNames = list(imagePath.glob(\"IMG_0000_*.tif\"))\n", + "panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "\n", + "if panelNames is not None:\n", + " panelCap = capture.Capture.from_filelist(panelNames)\n", + "else:\n", + " panelCap = None\n", + "\n", + "thecapture = capture.Capture.from_filelist(imageNames)\n", + "\n", + "# get camera model for future use\n", + "cam_model = thecapture.camera_model\n", + "# if this is a multicamera system like the RedEdge-MX Dual,\n", + "# we can combine the two serial numbers to help identify\n", + "# this camera system later.\n", + "if len(thecapture.camera_serials) > 1:\n", + " cam_serial = \"_\".join(thecapture.camera_serials)\n", + " print(cam_serial)\n", + "else:\n", + " cam_serial = thecapture.camera_serial\n", + "\n", + "print(\"Camera model:\", cam_model)\n", + "print(\"Bit depth:\", thecapture.bits_per_pixel)\n", + "print(\"Camera serial number:\", cam_serial)\n", + "print(\"Capture ID:\", thecapture.uuid)\n", + "\n", + "# determine if this sensor has a panchromatic band\n", + "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", + " panchroCam = True\n", + "else:\n", + " panchroCam = False\n", + " panSharpen = False\n", + "\n", + "if panelCap is not None:\n", + " if panelCap.panel_albedo() is not None:\n", + " panel_reflectance_by_band = panelCap.panel_albedo()\n", + " else:\n", + " panel_reflectance_by_band = [0.49] * len(\n", + " thecapture.eo_band_names()\n", + " ) # RedEdge band_index order\n", + " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)\n", + " irradiance_list = panelCap.panel_irradiance(panel_reflectance_by_band) + [\n", + " 0\n", + " ] # add to account for uncalibrated LWIR band, if applicable\n", + " img_type = \"reflectance\"\n", + " thecapture.plot_undistorted_reflectance(panel_irradiance)\n", + "else:\n", + " if thecapture.dls_present():\n", + " img_type = \"reflectance\"\n", + " irradiance_list = thecapture.dls_irradiance() + [0]\n", + " thecapture.plot_undistorted_reflectance(thecapture.dls_irradiance())\n", + " else:\n", + " img_type = \"radiance\"\n", + " thecapture.plot_undistorted_radiance()\n", + " irradiance_list = None" + ], + "execution_count": null, + "outputs": [], + "id": "7426068b" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Check for existing warp matrices \n", + "If we have already successfully aligned captures from this specific camera, we can typically save some time and use the alignment warp matrices for other captures from the same camera" + ], + "id": "8de5a774" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "from micasense.warp_io import load_warp_matrices\n", + "\n", + "if panchroCam:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", + "else:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", + "\n", + "if Path(\"./\" + warp_matrices_filename).is_file():\n", + " print(\"Found existing warp matrices for camera\", cam_serial)\n", + " loaded = load_warp_matrices(\n", + " warp_matrices_filename, as_projective=panchroCam\n", + " )\n", + " print(\"Warp matrices successfully loaded.\")\n", + "\n", + " if panchroCam:\n", + " warp_matrices_SIFT = loaded\n", + " else:\n", + " warp_matrices = loaded\n", + "else:\n", + " print(\"No existing warp matrices found. Create them later in the notebook.\")\n", + " warp_matrices_SIFT = False\n", + " warp_matrices = False" + ], + "execution_count": null, + "outputs": [], + "id": "35de1be0" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Unwarp and Align (OpenCV method for RedEdge3/M/MX/Dual and original Altum)\n", + "Alignment is a three-step process:\n", + "\n", + "Images are unwarped using the built-in lens calibration\n", + "A transformation is found to align each band to a common band\n", + "The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", + "We provide utilities to find the alignment transformations within a single capture. Our experience shows that once a good alignment 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 alignment 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." + ], + "id": "185627ff" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import cv2\n", + "import time\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", + "# We use a different alignment method for RedEdge-P and Altum-PT,\n", + "# so if the imagery is from this kind of camera, skip this step\n", + "if not panchroCam:\n", + " st = time.time()\n", + " # set to True if you'd like to ignore existing warp matrices and create new ones\n", + " regenerate = True\n", + " pyramid_levels = (\n", + " 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment\n", + " )\n", + " max_alignment_iterations = 10\n", + "\n", + " # match_index:\n", + " # for non-panchromatic cameras we want to use band 1, which is green\n", + " # the green band has zero rig relative offsets\n", + " # NOTE: These band numbers are zero-indexed\n", + " # special parameters for RedEdge-MX dual camera system\n", + " if len(thecapture.eo_band_names()) == 10:\n", + " print(\"is 10 band\")\n", + " pyramid_levels = 3\n", + " match_index = 4\n", + " max_alignment_iterations = 20\n", + " else:\n", + " match_index = 1\n", + "\n", + " warp_mode = (\n", + " cv2.MOTION_HOMOGRAPHY\n", + " ) # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", + "\n", + " print(\n", + " \"Aligning images. Depending on settings this can take from a few seconds to many minutes\"\n", + " )\n", + " # Can potentially increase max_iterations for better results, but longer runtimes\n", + " if warp_matrices and not regenerate:\n", + " print(\"Using existing warp matrices...\")\n", + " try:\n", + " irradiance = panel_irradiance + [0]\n", + " except NameError:\n", + " irradiance = None\n", + " else:\n", + " warp_matrices, alignment_pairs = imageutils.align_capture(\n", + " thecapture,\n", + " ref_index=match_index,\n", + " max_iterations=max_alignment_iterations,\n", + " warp_mode=warp_mode,\n", + " pyramid_levels=pyramid_levels,\n", + " )\n", + "\n", + " print(\"Finished Aligning\")\n", + " et = time.time()\n", + " elapsed_time = et - st\n", + " print(\"Alignment time:\", int(elapsed_time), \"seconds\")" + ], + "execution_count": null, + "outputs": [], + "id": "e154cdd1" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Crop Aligned Images and create aligned capture stack \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 take 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 to 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." + ], + "id": "c06bfdf3" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if not panchroCam:\n", + " cropped_dimensions, edges = imageutils.find_crop_bounds(\n", + " thecapture, warp_matrices, warp_mode=warp_mode, reference_band=match_index\n", + " )\n", + " print(\"Cropped dimensions:\", cropped_dimensions)\n", + " im_aligned = thecapture.create_aligned_capture(\n", + " warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type\n", + " )" + ], + "execution_count": null, + "outputs": [], + "id": "5fbfb346" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Alignment and pan-sharpening for RedEdge-P and Altum-PT\n", + "For older cameras we use OpenCV for capture alignment. For RedEdge-P and Altum-PT, we use SIFT (scale-invariant feature transform). We will use SIFT to create the warp matrices, then use these warp matrices to align the capture. See more here: \n", + "https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_sift.html\n", + "https://en.wikipedia.org/wiki/Scale-invariant_feature_transform\n", + "\n", + "For sensors with a panchromatic band (Altum-PT or RedEdge-P), we may wish to create a pan-sharpened stack. This example uses a linear interpolation method for pan-sharpening.\n", + "\n", + "The `radiometric_pan_sharpen` function takes in the SIFT warp matrices and outputs an upsampled image stack (all bands are changed to the resolution of the panchromatic sensor) and a pan-sharpened stack. \n", + "\n", + "Note: it is important to choose an initial alignment image that has a lot of straight, manmade features, such as roads or buildings. Trying to align a capture that is strictly vegetation will take a long time and may not be very good. This process can take a while, depending on how feature-rich the capture is. We have seen some Altum-PT captures align after 2 minutes, and some can take upwards of 30 minutes. Be patient! Once an alignment has been found, it can be used on other captures from the same camera, and the alignment process will be much faster. " + ], + "id": "c362623a" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "# from micasense.imageutils import brovey_pan_sharpen,radiometric_pan_sharpen\n", + "from skimage.transform import ProjectiveTransform\n", + "import time\n", + "\n", + "if panchroCam:\n", + " # set to True if you'd like to ignore existing warp matrices and create new ones\n", + " regenerate = True\n", + " st = time.time()\n", + " if not warp_matrices_SIFT or regenerate:\n", + " print(\"Generating new warp matrices...\")\n", + " warp_matrices_SIFT = thecapture.SIFT_align_capture(min_matches=10)\n", + "\n", + " sharpened_stack, upsampled = thecapture.radiometric_pan_sharpened_aligned_capture(\n", + " warp_matrices=warp_matrices_SIFT,\n", + " irradiance_list=irradiance_list,\n", + " img_type=img_type,\n", + " )\n", + "\n", + " # we can also use the Rig Relatives from the image metadata to do a quick, rudimentary alignment\n", + " # warp_matrices0=thecapture.get_warp_matrices(ref_index=5)\n", + " # sharpened_stack,upsampled = radiometric_pan_sharpen(thecapture,warp_matrices=warp_matrices0)\n", + "\n", + " print(\"Pansharpened shape:\", sharpened_stack.shape)\n", + " print(\"Upsampled shape:\", upsampled.shape)\n", + " # re-assign to im_aligned to match rest of code\n", + " im_aligned = upsampled\n", + " et = time.time()\n", + " elapsed_time = et - st\n", + " print(\"Alignment and pan-sharpening time:\", int(elapsed_time), \"seconds\")" + ], + "execution_count": null, + "outputs": [], + "id": "12c2a278" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Save warp matrices\n", + "Once an alignment for your camera has been found, it can be saved to a file for later use with this notebook. It can also be used on the Batch Alignment notebook for aligning all of the captures from an entire flight." + ], + "id": "3e9d51b1" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "from micasense.warp_io import save_warp_matrices\n", + "\n", + "if panchroCam:\n", + " working_wm = warp_matrices_SIFT\n", + "else:\n", + " working_wm = warp_matrices\n", + "if not Path(\"./\" + warp_matrices_filename).is_file() or regenerate:\n", + " save_warp_matrices(warp_matrices_filename, working_wm)\n", + " print(\"Saved to\", Path(\"./\" + warp_matrices_filename).resolve())\n", + "else:\n", + " print(\"Matrices already exist at\", Path(\"./\" + warp_matrices_filename).resolve())" + ], + "execution_count": null, + "outputs": [], + "id": "84855db4" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Multispectral band histogram \n", + "We can compare the radiance between multispectral bands, and between upsampled/pansharpened in the case of RedEdge-P and Altum-PT" + ], + "id": "47d3fb77" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "theColors = {\n", + " \"Blue\": \"blue\",\n", + " \"Green\": \"green\",\n", + " \"Red\": \"red\",\n", + " \"Red edge\": \"maroon\",\n", + " \"NIR\": \"purple\",\n", + " \"Panchro\": \"yellow\",\n", + " \"PanchroB\": \"orange\",\n", + " \"Red edge-740\": \"salmon\",\n", + " \"Red Edge\": \"maroon\",\n", + " \"Blue-444\": \"aqua\",\n", + " \"Green-531\": \"lime\",\n", + " \"Red-650\": \"lightcoral\",\n", + " \"Red edge-705\": \"brown\",\n", + "}\n", + "\n", + "eo_count = len(thecapture.eo_indices())\n", + "multispec_min = np.min(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 0.01))\n", + "multispec_max = np.max(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 99.99))\n", + "\n", + "theRange = (multispec_min, multispec_max)\n", + "\n", + "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + "for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", + " axis.hist(\n", + " im_aligned[:, :, x].ravel(),\n", + " bins=512,\n", + " range=theRange,\n", + " histtype=\"step\",\n", + " label=y,\n", + " color=theColors[y],\n", + " linewidth=1.5,\n", + " )\n", + "plt.title(\"Multispectral histogram (radiance)\")\n", + "axis.legend()\n", + "plt.show()\n", + "\n", + "if panchroCam:\n", + " eo_count = len(thecapture.eo_indices())\n", + " multispec_min = np.min(\n", + " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 0.01)\n", + " )\n", + " multispec_max = np.max(\n", + " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 99.99)\n", + " )\n", + "\n", + " theRange = (multispec_min, multispec_max)\n", + "\n", + " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + " for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", + " axis.hist(\n", + " sharpened_stack[:, :, x].ravel(),\n", + " bins=512,\n", + " range=theRange,\n", + " histtype=\"step\",\n", + " label=y,\n", + " color=theColors[y],\n", + " linewidth=1.5,\n", + " )\n", + " plt.title(\"Pan-sharpened multispectral histogram (radiance)\")\n", + " axis.legend()\n", + " plt.show()" + ], + "execution_count": null, + "outputs": [], + "id": "8e036b45" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualize Aligned Images\n", + "Once the transformation has been found, it can be verified by compositing 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 external software such as QGIS. Useful composites are a naturally colored RGB as well as color infrared, or CIR." + ], + "id": "3ccb819b" + }, + { + "cell_type": "code", + "metadata": {}, + "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 = [\n", + " thecapture.band_names_lower().index(\"red\"),\n", + " thecapture.band_names_lower().index(\"green\"),\n", + " thecapture.band_names_lower().index(\"blue\"),\n", + "]\n", + "cir_band_indices = [\n", + " thecapture.band_names_lower().index(\"nir\"),\n", + " thecapture.band_names_lower().index(\"red\"),\n", + " thecapture.band_names_lower().index(\"green\"),\n", + "]\n", + "\n", + "# Create normalized stacks for viewing\n", + "im_display = np.zeros(\n", + " (im_aligned.shape[0], im_aligned.shape[1], im_aligned.shape[2]), dtype=float\n", + ")\n", + "im_min = np.percentile(\n", + " im_aligned[:, :, rgb_band_indices].flatten(), 0.5\n", + ") # modify these percentiles to adjust contrast\n", + "im_max = np.percentile(\n", + " im_aligned[:, :, rgb_band_indices].flatten(), 99.5\n", + ") # for many images, 0.5 and 99.5 are good values\n", + "\n", + "if panchroCam:\n", + " im_display_sharp = np.zeros(\n", + " (sharpened_stack.shape[0], sharpened_stack.shape[1], sharpened_stack.shape[2]),\n", + " dtype=float,\n", + " )\n", + " im_min_sharp = np.percentile(\n", + " sharpened_stack[:, :, rgb_band_indices].flatten(), 0.5\n", + " ) # modify these percentiles to adjust contrast\n", + " im_max_sharp = np.percentile(\n", + " sharpened_stack[:, :, rgb_band_indices].flatten(), 99.5\n", + " ) # for many images, 0.5 and 99.5 are good values\n", + "\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", + " if panchroCam:\n", + " im_display_sharp[:, :, i] = imageutils.normalize(\n", + " sharpened_stack[:, :, i], im_min_sharp, im_max_sharp\n", + " )\n", + "\n", + "rgb = im_display[:, :, rgb_band_indices]\n", + "\n", + "if panchroCam:\n", + " rgb_sharp = im_display_sharp[:, :, rgb_band_indices]\n", + "\n", + "nir_band = thecapture.band_names_lower().index(\"nir\")\n", + "red_band = thecapture.band_names_lower().index(\"red\")\n", + "\n", + "ndvi = (im_aligned[:, :, nir_band] - im_aligned[:, :, red_band]) / (\n", + " im_aligned[:, :, nir_band] + im_aligned[:, :, red_band]\n", + ")\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", + "if panchroCam:\n", + " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)\n", + "else:\n", + " fig, ax1 = plt.subplots(1, 1, figsize=figsize)\n", + "ax1.set_title(\"Red-Green-Blue Composite\")\n", + "ax1.imshow(rgb)\n", + "if panchroCam:\n", + " ax2.set_title(\"Red-Green-Blue Composite (pan-sharpened)\")\n", + " ax2.imshow(rgb_sharp)\n", + "\n", + "fig, (ax3, ax4) = plt.subplots(1, 2, figsize=figsize)\n", + "ax3.set_title(\"NDVI\")\n", + "ax3.imshow(ndvi)\n", + "ax4.set_title(\"Color Infrared (CIR) Composite\")\n", + "ax4.imshow(cir)\n", + "\n", + "# set custom lims if you want to zoom in to image to see more detail\n", + "# this is useful for comparing upsampled and pan-sharpened stacks\n", + "# custom_xlim=(1500,2000)\n", + "# custom_ylim=(2000,1500)\n", + "# plt.setp([ax1,ax2,ax3,ax4], xlim=custom_xlim, ylim=custom_ylim)\n", + "\n", + "plt.show()" + ], + "execution_count": null, + "outputs": [], + "id": "8626c90b" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Image Enhancement\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." + ], + "id": "17f71519" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if panchroCam:\n", + " rgb = rgb_sharp\n", + "# 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()" + ], + "execution_count": null, + "outputs": [], + "id": "1185fd16" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Stack Export\n", + "We can easily export the stacks into an image 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. The stacks include geographic information. \n", + "\n", + "If you prefer, you may set `sort_by_wavelength` to `True` in the `save_capture_as_stack` function.\n", + "\n", + "Unless otherwise specified, this will save in your working folder, that is your `imageprocessing` directory." + ], + "id": "6d1f6189" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "# set output name to unique capture ID, e.g. FWoNSvgDNBX63Xv378qs\n", + "outputName = thecapture.uuid\n", + "\n", + "st = time.time()\n", + "if panchroCam:\n", + " # in this example, we can export both a pan-sharpened stack and an upsampled stack\n", + " # so you can compare them in GIS. In practice, you would typically only output the pansharpened stack\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-pansharpened.tif\", sort_by_wavelength=True, pansharpen=True\n", + " )\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-upsampled.tif\", sort_by_wavelength=True, pansharpen=False\n", + " )\n", + "else:\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-noPanels.tif\", sort_by_wavelength=True\n", + " )\n", + "\n", + "et = time.time()\n", + "elapsed_time = et - st\n", + "print(\"Time to save stacks:\", int(elapsed_time), \"seconds.\")" + ], + "execution_count": null, + "outputs": [], + "id": "a2d716be" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# NDVI Computation\n", + "For raw index computation on single images, the `numpy` package provides a simple way to do math and simple visualization on images. Below, we compute and visualize an image histogram, and then use that to pick a color map range for visualizing the NDVI of an image.\n", + "\n", + "## Plant Classification\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." + ], + "id": "f4d9cd64" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "nir_band = thecapture.band_names_lower().index(\"nir\")\n", + "red_band = thecapture.band_names_lower().index(\"red\")\n", + "\n", + "thelayer = im_aligned\n", + "if panchroCam:\n", + " thelayer = sharpened_stack\n", + "np.seterr(\n", + " divide=\"ignore\", invalid=\"ignore\"\n", + ") # 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 = (thelayer[:, :, nir_band] - thelayer[:, :, red_band]) / (\n", + " thelayer[:, :, nir_band] + thelayer[:, :, red_band]\n", + ")\n", + "print(\"Image type:\", img_type)\n", + "\n", + "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", + "# this does not seem to work on panchro stacks\n", + "if img_type == \"reflectance\":\n", + " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < 0.20, ndvi)\n", + "elif img_type == \"radiance\":\n", + " lower_pct_radiance = np.percentile(thelayer[:, :, nir_band], 10.0)\n", + " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < lower_pct_radiance, ndvi)\n", + "# Compute and display a histogram\n", + "# ndvi_hist_min = np.min(ndvi)\n", + "# ndvi_hist_max = np.max(ndvi)\n", + "ndvi_hist_min = np.min(np.percentile(ndvi, 0.5))\n", + "ndvi_hist_max = np.max(np.percentile(ndvi, 99.5))\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(\n", + " ndvi.flatten(), 99.5\n", + ") # 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(\n", + " gamma_corr_rgb,\n", + " masked_ndvi,\n", + " figsize=(14, 7),\n", + " title=\"NDVI filtered to only plants over RGB base layer\",\n", + " vmin=min_display_ndvi,\n", + " vmax=max_display_ndvi,\n", + ")\n", + "fig.savefig(thecapture.uuid + \"_ndvi_over_rgb.png\")" + ], + "execution_count": null, + "outputs": [], + "id": "4a572a55" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# NDRE Computation\n", + "In the same manner, we can compute, filter, and display another index useful for MicaSense cameras, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health." + ], + "id": "d8b4cf19" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n", + "rededge_band = thecapture.band_names_lower().index(\"red edge\")\n", + "ndre = (thelayer[:, :, nir_band] - thelayer[:, :, rededge_band]) / (\n", + " thelayer[:, :, nir_band] + thelayer[:, :, rededge_band]\n", + ")\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(np.percentile(masked_ndre, 0.5))\n", + "ndre_hist_max = np.max(np.percentile(masked_ndre, 99.5))\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(\n", + " gamma_corr_rgb,\n", + " masked_ndre,\n", + " figsize=(14, 7),\n", + " title=\"NDRE filtered to only plants over RGB base layer\",\n", + " vmin=min_display_ndre,\n", + " vmax=max_display_ndre,\n", + ")\n", + "fig.savefig(thecapture.uuid + \"_ndre_over_rgb.png\")" + ], + "execution_count": null, + "outputs": [], + "id": "b3fe7a83" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Thermal Imagery\n", + "If our image is from an Altum or Altum-PT 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" + ], + "id": "898cc28a" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if len(thecapture.lw_indices()) > 0:\n", + " lwir_band = thecapture.band_names_lower().index(\"lwir\")\n", + "\n", + " # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n", + " if panchroCam:\n", + " masked_thermal = sharpened_stack[:, :, lwir_band]\n", + " else:\n", + " masked_thermal = im_aligned[:, :, lwir_band]\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", + " # Compute a histogram\n", + " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + " axis.hist(\n", + " masked_thermal.ravel(),\n", + " bins=512,\n", + " range=(\n", + " np.min(np.percentile(masked_thermal, 1)),\n", + " np.max(np.percentile(masked_thermal, 99)),\n", + " ),\n", + " )\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(\n", + " gamma_corr_rgb,\n", + " masked_thermal,\n", + " figsize=(14, 7),\n", + " title=\"Temperature over True Color\",\n", + " vmin=min_display_therm,\n", + " 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=0.4,\n", + " contour_fmt=\"%.0fC\",\n", + " )\n", + " fig.savefig(thecapture.uuid + \"_thermal_over_rgb.png\")" + ], + "execution_count": null, + "outputs": [], + "id": "176e2c34" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Red vs NIR Reflectance\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." + ], + "id": "bc9b29ec" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "x_band = red_band\n", + "y_band = nir_band\n", + "x_max = np.max(np.percentile(im_aligned[:, :, x_band], 99.99))\n", + "y_max = np.max(np.percentile(im_aligned[:, :, y_band], 99.99))\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "plt.hexbin(\n", + " im_aligned[:, :, x_band],\n", + " im_aligned[:, :, y_band],\n", + " gridsize=640,\n", + " bins=\"log\",\n", + " extent=(0, x_max, 0, y_max),\n", + ")\n", + "ax = fig.gca()\n", + "ax.set_xlim([0, x_max])\n", + "ax.set_ylim([0, y_max])\n", + "plt.xlabel(\"{} Reflectance\".format(thecapture.band_names()[x_band]))\n", + "plt.ylabel(\"{} Reflectance\".format(thecapture.band_names()[y_band]))\n", + "plt.show()\n", + "\n", + "if panchroCam:\n", + " x_band = red_band\n", + " y_band = nir_band\n", + " x_max = np.max(np.percentile(sharpened_stack[:, :, x_band], 99.99))\n", + " y_max = np.max(np.percentile(sharpened_stack[:, :, y_band], 99.99))\n", + "\n", + " fig = plt.figure(figsize=(12, 12))\n", + " plt.hexbin(\n", + " sharpened_stack[:, :, x_band],\n", + " sharpened_stack[:, :, y_band],\n", + " gridsize=640,\n", + " bins=\"log\",\n", + " extent=(0, x_max, 0, y_max),\n", + " )\n", + " ax = fig.gca()\n", + " ax.set_xlim([0, x_max])\n", + " ax.set_ylim([0, y_max])\n", + " plt.xlabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[x_band]))\n", + " plt.ylabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[y_band]))\n", + " plt.show()" + ], + "execution_count": null, + "outputs": [], + "id": "c0250872" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Print warp_matrices for usage elsewhere, such as Batch Processing\n", + "Lastly, 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." + ], + "id": "4ec2d106" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if panchroCam:\n", + " print(warp_matrices_SIFT)\n", + "else:\n", + " print(warp_matrices)" + ], + "execution_count": null, + "outputs": [], + "id": "c7cf8f48" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.12" + } }, - { - "cell_type": "code", - "execution_count": null, - "id": "b09cd3b1", - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7426068b", - "metadata": {}, - "outputs": [], - "source": [ - "import micasense.capture as capture\n", - "\n", - "%matplotlib inline\n", - "from pathlib import Path\n", - "import matplotlib.pyplot as plt\n", - "\n", - "plt.rcParams[\"figure.facecolor\"] = \"w\"\n", - "\n", - "panelNames = None\n", - "\n", - "# set your image paths here. See more here: https://docs.python.org/3/library/pathlib.html\n", - "# if using Windows, you need to an an \"r\" to the path like this: Path(r\"C:\\Files\")\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-MX-DUAL\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0431_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-MX\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0020_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0001_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/ALTUM\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0021_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-P\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0011_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "imagePath = Path(\"./data/ALTUM-PT\")\n", - "\n", - "# these will return lists of image paths as strings\n", - "imageNames = list(imagePath.glob(\"IMG_0010_*.tif\"))\n", - "imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "panelNames = list(imagePath.glob(\"IMG_0000_*.tif\"))\n", - "panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "\n", - "if panelNames is not None:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "else:\n", - " panelCap = None\n", - "\n", - "thecapture = capture.Capture.from_filelist(imageNames)\n", - "\n", - "# get camera model for future use\n", - "cam_model = thecapture.camera_model\n", - "# if this is a multicamera system like the RedEdge-MX Dual,\n", - "# we can combine the two serial numbers to help identify\n", - "# this camera system later.\n", - "if len(thecapture.camera_serials) > 1:\n", - " cam_serial = \"_\".join(thecapture.camera_serials)\n", - " print(cam_serial)\n", - "else:\n", - " cam_serial = thecapture.camera_serial\n", - "\n", - "print(\"Camera model:\", cam_model)\n", - "print(\"Bit depth:\", thecapture.bits_per_pixel)\n", - "print(\"Camera serial number:\", cam_serial)\n", - "print(\"Capture ID:\", thecapture.uuid)\n", - "\n", - "# determine if this sensor has a panchromatic band\n", - "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", - " panchroCam = True\n", - "else:\n", - " panchroCam = False\n", - " panSharpen = False\n", - "\n", - "if panelCap is not None:\n", - " if panelCap.panel_albedo() is not None:\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " panel_reflectance_by_band = [0.49] * len(\n", - " thecapture.eo_band_names()\n", - " ) # RedEdge band_index order\n", - " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)\n", - " irradiance_list = panelCap.panel_irradiance(panel_reflectance_by_band) + [\n", - " 0\n", - " ] # add to account for uncalibrated LWIR band, if applicable\n", - " img_type = \"reflectance\"\n", - " thecapture.plot_undistorted_reflectance(panel_irradiance)\n", - "else:\n", - " if thecapture.dls_present():\n", - " img_type = \"reflectance\"\n", - " irradiance_list = thecapture.dls_irradiance() + [0]\n", - " thecapture.plot_undistorted_reflectance(thecapture.dls_irradiance())\n", - " else:\n", - " img_type = \"radiance\"\n", - " thecapture.plot_undistorted_radiance()\n", - " irradiance_list = None" - ] - }, - { - "cell_type": "markdown", - "id": "8de5a774", - "metadata": {}, - "source": [ - "# Check for existing warp matrices \n", - "If we have already successfully aligned captures from this specific camera, we can typically save some time and use the alignment warp matrices for other captures from the same camera" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "35de1be0", - "metadata": {}, - "outputs": [], - "source": [ - "from skimage.transform import ProjectiveTransform\n", - "import numpy as np\n", - "\n", - "if panchroCam:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", - "else:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", - "\n", - "if Path(\"./\" + warp_matrices_filename).is_file():\n", - " print(\"Found existing warp matrices for camera\", cam_serial)\n", - " load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True)\n", - " loaded_warp_matrices = []\n", - " for matrix in load_warp_matrices:\n", - " if panchroCam:\n", - " transform = ProjectiveTransform(matrix=matrix.astype(\"float64\"))\n", - " loaded_warp_matrices.append(transform)\n", - " else:\n", - " loaded_warp_matrices.append(matrix.astype(\"float32\"))\n", - " print(\"Warp matrices successfully loaded.\")\n", - "\n", - " if panchroCam:\n", - " warp_matrices_SIFT = loaded_warp_matrices\n", - " else:\n", - " warp_matrices = loaded_warp_matrices\n", - "else:\n", - " print(\"No existing warp matrices found. Create them later in the notebook.\")\n", - " warp_matrices_SIFT = False\n", - " warp_matrices = False" - ] - }, - { - "cell_type": "markdown", - "id": "185627ff", - "metadata": {}, - "source": [ - "# Unwarp and Align (OpenCV method for RedEdge3/M/MX/Dual and original Altum)\n", - "Alignment is a three-step process:\n", - "\n", - "Images are unwarped using the built-in lens calibration\n", - "A transformation is found to align each band to a common band\n", - "The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", - "We provide utilities to find the alignment transformations within a single capture. Our experience shows that once a good alignment 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 alignment 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, - "id": "e154cdd1", - "metadata": {}, - "outputs": [], - "source": [ - "import cv2\n", - "import time\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", - "# We use a different alignment method for RedEdge-P and Altum-PT,\n", - "# so if the imagery is from this kind of camera, skip this step\n", - "if not panchroCam:\n", - " st = time.time()\n", - " # set to True if you'd like to ignore existing warp matrices and create new ones\n", - " regenerate = True\n", - " pyramid_levels = (\n", - " 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment\n", - " )\n", - " max_alignment_iterations = 10\n", - "\n", - " # match_index:\n", - " # for non-panchromatic cameras we want to use band 1, which is green\n", - " # the green band has zero rig relative offsets\n", - " # NOTE: These band numbers are zero-indexed\n", - " # special parameters for RedEdge-MX dual camera system\n", - " if len(thecapture.eo_band_names()) == 10:\n", - " print(\"is 10 band\")\n", - " pyramid_levels = 3\n", - " match_index = 4\n", - " max_alignment_iterations = 20\n", - " else:\n", - " match_index = 1\n", - "\n", - " warp_mode = (\n", - " cv2.MOTION_HOMOGRAPHY\n", - " ) # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", - "\n", - " print(\n", - " \"Aligning images. Depending on settings this can take from a few seconds to many minutes\"\n", - " )\n", - " # Can potentially increase max_iterations for better results, but longer runtimes\n", - " if warp_matrices and not regenerate:\n", - " print(\"Using existing warp matrices...\")\n", - " try:\n", - " irradiance = panel_irradiance + [0]\n", - " except NameError:\n", - " irradiance = None\n", - " else:\n", - " warp_matrices, alignment_pairs = imageutils.align_capture(\n", - " thecapture,\n", - " ref_index=match_index,\n", - " max_iterations=max_alignment_iterations,\n", - " warp_mode=warp_mode,\n", - " pyramid_levels=pyramid_levels,\n", - " )\n", - "\n", - " print(\"Finished Aligning\")\n", - " et = time.time()\n", - " elapsed_time = et - st\n", - " print(\"Alignment time:\", int(elapsed_time), \"seconds\")" - ] - }, - { - "cell_type": "markdown", - "id": "c06bfdf3", - "metadata": {}, - "source": [ - "# Crop Aligned Images and create aligned capture stack \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 take 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 to 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, - "id": "5fbfb346", - "metadata": {}, - "outputs": [], - "source": [ - "if not panchroCam:\n", - " cropped_dimensions, edges = imageutils.find_crop_bounds(\n", - " thecapture, warp_matrices, warp_mode=warp_mode, reference_band=match_index\n", - " )\n", - " print(\"Cropped dimensions:\", cropped_dimensions)\n", - " im_aligned = thecapture.create_aligned_capture(\n", - " warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type\n", - " )" - ] - }, - { - "cell_type": "markdown", - "id": "c362623a", - "metadata": {}, - "source": [ - "# Alignment and pan-sharpening for RedEdge-P and Altum-PT\n", - "For older cameras we use OpenCV for capture alignment. For RedEdge-P and Altum-PT, we use SIFT (scale-invariant feature transform). We will use SIFT to create the warp matrices, then use these warp matrices to align the capture. See more here: \n", - "https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_sift.html\n", - "https://en.wikipedia.org/wiki/Scale-invariant_feature_transform\n", - "\n", - "For sensors with a panchromatic band (Altum-PT or RedEdge-P), we may wish to create a pan-sharpened stack. This example uses a linear interpolation method for pan-sharpening.\n", - "\n", - "The `radiometric_pan_sharpen` function takes in the SIFT warp matrices and outputs an upsampled image stack (all bands are changed to the resolution of the panchromatic sensor) and a pan-sharpened stack. \n", - "\n", - "Note: it is important to choose an initial alignment image that has a lot of straight, manmade features, such as roads or buildings. Trying to align a capture that is strictly vegetation will take a long time and may not be very good. This process can take a while, depending on how feature-rich the capture is. We have seen some Altum-PT captures align after 2 minutes, and some can take upwards of 30 minutes. Be patient! Once an alignment has been found, it can be used on other captures from the same camera, and the alignment process will be much faster. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "12c2a278", - "metadata": {}, - "outputs": [], - "source": [ - "# from micasense.imageutils import brovey_pan_sharpen,radiometric_pan_sharpen\n", - "from skimage.transform import ProjectiveTransform\n", - "import time\n", - "\n", - "if panchroCam:\n", - " # set to True if you'd like to ignore existing warp matrices and create new ones\n", - " regenerate = True\n", - " st = time.time()\n", - " if not warp_matrices_SIFT or regenerate:\n", - " print(\"Generating new warp matrices...\")\n", - " warp_matrices_SIFT = thecapture.SIFT_align_capture(min_matches=10)\n", - "\n", - " sharpened_stack, upsampled = thecapture.radiometric_pan_sharpened_aligned_capture(\n", - " warp_matrices=warp_matrices_SIFT,\n", - " irradiance_list=irradiance_list,\n", - " img_type=img_type,\n", - " )\n", - "\n", - " # we can also use the Rig Relatives from the image metadata to do a quick, rudimentary alignment\n", - " # warp_matrices0=thecapture.get_warp_matrices(ref_index=5)\n", - " # sharpened_stack,upsampled = radiometric_pan_sharpen(thecapture,warp_matrices=warp_matrices0)\n", - "\n", - " print(\"Pansharpened shape:\", sharpened_stack.shape)\n", - " print(\"Upsampled shape:\", upsampled.shape)\n", - " # re-assign to im_aligned to match rest of code\n", - " im_aligned = upsampled\n", - " et = time.time()\n", - " elapsed_time = et - st\n", - " print(\"Alignment and pan-sharpening time:\", int(elapsed_time), \"seconds\")" - ] - }, - { - "cell_type": "markdown", - "id": "3e9d51b1", - "metadata": {}, - "source": [ - "# Save warp matrices\n", - "Once an alignment for your camera has been found, it can be saved to a file for later use with this notebook. It can also be used on the Batch Alignment notebook for aligning all of the captures from an entire flight." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "84855db4", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy\n", - "import skimage\n", - "from skimage.transform import ProjectiveTransform\n", - "\n", - "if panchroCam:\n", - " working_wm = warp_matrices_SIFT\n", - "else:\n", - " working_wm = warp_matrices\n", - "if not Path(\"./\" + warp_matrices_filename).is_file() or regenerate:\n", - " temp_matrices = []\n", - " for x in working_wm:\n", - " if isinstance(x, numpy.ndarray):\n", - " temp_matrices.append(x)\n", - " if isinstance(x, skimage.transform._geometric.ProjectiveTransform):\n", - " temp_matrices.append(x.params)\n", - " np.save(\n", - " warp_matrices_filename, np.array(temp_matrices, dtype=object), allow_pickle=True\n", - " )\n", - " print(\"Saved to\", Path(\"./\" + warp_matrices_filename).resolve())\n", - "else:\n", - " print(\"Matrices already exist at\", Path(\"./\" + warp_matrices_filename).resolve())" - ] - }, - { - "cell_type": "markdown", - "id": "47d3fb77", - "metadata": {}, - "source": [ - "# Multispectral band histogram \n", - "We can compare the radiance between multispectral bands, and between upsampled/pansharpened in the case of RedEdge-P and Altum-PT" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8e036b45", - "metadata": {}, - "outputs": [], - "source": [ - "theColors = {\n", - " \"Blue\": \"blue\",\n", - " \"Green\": \"green\",\n", - " \"Red\": \"red\",\n", - " \"Red edge\": \"maroon\",\n", - " \"NIR\": \"purple\",\n", - " \"Panchro\": \"yellow\",\n", - " \"PanchroB\": \"orange\",\n", - " \"Red edge-740\": \"salmon\",\n", - " \"Red Edge\": \"maroon\",\n", - " \"Blue-444\": \"aqua\",\n", - " \"Green-531\": \"lime\",\n", - " \"Red-650\": \"lightcoral\",\n", - " \"Red edge-705\": \"brown\",\n", - "}\n", - "\n", - "eo_count = len(thecapture.eo_indices())\n", - "multispec_min = np.min(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 0.01))\n", - "multispec_max = np.max(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 99.99))\n", - "\n", - "theRange = (multispec_min, multispec_max)\n", - "\n", - "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - "for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", - " axis.hist(\n", - " im_aligned[:, :, x].ravel(),\n", - " bins=512,\n", - " range=theRange,\n", - " histtype=\"step\",\n", - " label=y,\n", - " color=theColors[y],\n", - " linewidth=1.5,\n", - " )\n", - "plt.title(\"Multispectral histogram (radiance)\")\n", - "axis.legend()\n", - "plt.show()\n", - "\n", - "if panchroCam:\n", - " eo_count = len(thecapture.eo_indices())\n", - " multispec_min = np.min(\n", - " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 0.01)\n", - " )\n", - " multispec_max = np.max(\n", - " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 99.99)\n", - " )\n", - "\n", - " theRange = (multispec_min, multispec_max)\n", - "\n", - " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - " for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", - " axis.hist(\n", - " sharpened_stack[:, :, x].ravel(),\n", - " bins=512,\n", - " range=theRange,\n", - " histtype=\"step\",\n", - " label=y,\n", - " color=theColors[y],\n", - " linewidth=1.5,\n", - " )\n", - " plt.title(\"Pan-sharpened multispectral histogram (radiance)\")\n", - " axis.legend()\n", - " plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "3ccb819b", - "metadata": {}, - "source": [ - "# Visualize Aligned Images\n", - "Once the transformation has been found, it can be verified by compositing 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 external software such as QGIS. Useful composites are a naturally colored RGB as well as color infrared, or CIR." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8626c90b", - "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 = [\n", - " thecapture.band_names_lower().index(\"red\"),\n", - " thecapture.band_names_lower().index(\"green\"),\n", - " thecapture.band_names_lower().index(\"blue\"),\n", - "]\n", - "cir_band_indices = [\n", - " thecapture.band_names_lower().index(\"nir\"),\n", - " thecapture.band_names_lower().index(\"red\"),\n", - " thecapture.band_names_lower().index(\"green\"),\n", - "]\n", - "\n", - "# Create normalized stacks for viewing\n", - "im_display = np.zeros(\n", - " (im_aligned.shape[0], im_aligned.shape[1], im_aligned.shape[2]), dtype=float\n", - ")\n", - "im_min = np.percentile(\n", - " im_aligned[:, :, rgb_band_indices].flatten(), 0.5\n", - ") # modify these percentiles to adjust contrast\n", - "im_max = np.percentile(\n", - " im_aligned[:, :, rgb_band_indices].flatten(), 99.5\n", - ") # for many images, 0.5 and 99.5 are good values\n", - "\n", - "if panchroCam:\n", - " im_display_sharp = np.zeros(\n", - " (sharpened_stack.shape[0], sharpened_stack.shape[1], sharpened_stack.shape[2]),\n", - " dtype=float,\n", - " )\n", - " im_min_sharp = np.percentile(\n", - " sharpened_stack[:, :, rgb_band_indices].flatten(), 0.5\n", - " ) # modify these percentiles to adjust contrast\n", - " im_max_sharp = np.percentile(\n", - " sharpened_stack[:, :, rgb_band_indices].flatten(), 99.5\n", - " ) # for many images, 0.5 and 99.5 are good values\n", - "\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", - " if panchroCam:\n", - " im_display_sharp[:, :, i] = imageutils.normalize(\n", - " sharpened_stack[:, :, i], im_min_sharp, im_max_sharp\n", - " )\n", - "\n", - "rgb = im_display[:, :, rgb_band_indices]\n", - "\n", - "if panchroCam:\n", - " rgb_sharp = im_display_sharp[:, :, rgb_band_indices]\n", - "\n", - "nir_band = thecapture.band_names_lower().index(\"nir\")\n", - "red_band = thecapture.band_names_lower().index(\"red\")\n", - "\n", - "ndvi = (im_aligned[:, :, nir_band] - im_aligned[:, :, red_band]) / (\n", - " im_aligned[:, :, nir_band] + im_aligned[:, :, red_band]\n", - ")\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", - "if panchroCam:\n", - " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)\n", - "else:\n", - " fig, ax1 = plt.subplots(1, 1, figsize=figsize)\n", - "ax1.set_title(\"Red-Green-Blue Composite\")\n", - "ax1.imshow(rgb)\n", - "if panchroCam:\n", - " ax2.set_title(\"Red-Green-Blue Composite (pan-sharpened)\")\n", - " ax2.imshow(rgb_sharp)\n", - "\n", - "fig, (ax3, ax4) = plt.subplots(1, 2, figsize=figsize)\n", - "ax3.set_title(\"NDVI\")\n", - "ax3.imshow(ndvi)\n", - "ax4.set_title(\"Color Infrared (CIR) Composite\")\n", - "ax4.imshow(cir)\n", - "\n", - "# set custom lims if you want to zoom in to image to see more detail\n", - "# this is useful for comparing upsampled and pan-sharpened stacks\n", - "# custom_xlim=(1500,2000)\n", - "# custom_ylim=(2000,1500)\n", - "# plt.setp([ax1,ax2,ax3,ax4], xlim=custom_xlim, ylim=custom_ylim)\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "17f71519", - "metadata": {}, - "source": [ - "# Image Enhancement\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, - "id": "1185fd16", - "metadata": {}, - "outputs": [], - "source": [ - "if panchroCam:\n", - " rgb = rgb_sharp\n", - "# 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", - "id": "6d1f6189", - "metadata": {}, - "source": [ - "# Stack Export\n", - "We can easily export the stacks into an image 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. The stacks include geographic information. \n", - "\n", - "If you prefer, you may set `sort_by_wavelength` to `True` in the `save_capture_as_stack` function.\n", - "\n", - "Unless otherwise specified, this will save in your working folder, that is your `imageprocessing` directory." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a2d716be", - "metadata": {}, - "outputs": [], - "source": [ - "# set output name to unique capture ID, e.g. FWoNSvgDNBX63Xv378qs\n", - "outputName = thecapture.uuid\n", - "\n", - "st = time.time()\n", - "if panchroCam:\n", - " # in this example, we can export both a pan-sharpened stack and an upsampled stack\n", - " # so you can compare them in GIS. In practice, you would typically only output the pansharpened stack\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-pansharpened.tif\", sort_by_wavelength=True, pansharpen=True\n", - " )\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-upsampled.tif\", sort_by_wavelength=True, pansharpen=False\n", - " )\n", - "else:\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-noPanels.tif\", sort_by_wavelength=True\n", - " )\n", - "\n", - "et = time.time()\n", - "elapsed_time = et - st\n", - "print(\"Time to save stacks:\", int(elapsed_time), \"seconds.\")" - ] - }, - { - "cell_type": "markdown", - "id": "f4d9cd64", - "metadata": {}, - "source": [ - "# NDVI Computation\n", - "For raw index computation on single images, the `numpy` package provides a simple way to do math and simple visualization on images. Below, we compute and visualize an image histogram, and then use that to pick a color map range for visualizing the NDVI of an image.\n", - "\n", - "## Plant Classification\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, - "id": "4a572a55", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "nir_band = thecapture.band_names_lower().index(\"nir\")\n", - "red_band = thecapture.band_names_lower().index(\"red\")\n", - "\n", - "thelayer = im_aligned\n", - "if panchroCam:\n", - " thelayer = sharpened_stack\n", - "np.seterr(\n", - " divide=\"ignore\", invalid=\"ignore\"\n", - ") # 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 = (thelayer[:, :, nir_band] - thelayer[:, :, red_band]) / (\n", - " thelayer[:, :, nir_band] + thelayer[:, :, red_band]\n", - ")\n", - "print(\"Image type:\", img_type)\n", - "\n", - "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", - "# this does not seem to work on panchro stacks\n", - "if img_type == \"reflectance\":\n", - " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < 0.20, ndvi)\n", - "elif img_type == \"radiance\":\n", - " lower_pct_radiance = np.percentile(thelayer[:, :, nir_band], 10.0)\n", - " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < lower_pct_radiance, ndvi)\n", - "# Compute and display a histogram\n", - "# ndvi_hist_min = np.min(ndvi)\n", - "# ndvi_hist_max = np.max(ndvi)\n", - "ndvi_hist_min = np.min(np.percentile(ndvi, 0.5))\n", - "ndvi_hist_max = np.max(np.percentile(ndvi, 99.5))\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(\n", - " ndvi.flatten(), 99.5\n", - ") # 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(\n", - " gamma_corr_rgb,\n", - " masked_ndvi,\n", - " figsize=(14, 7),\n", - " title=\"NDVI filtered to only plants over RGB base layer\",\n", - " vmin=min_display_ndvi,\n", - " vmax=max_display_ndvi,\n", - ")\n", - "fig.savefig(thecapture.uuid + \"_ndvi_over_rgb.png\")" - ] - }, - { - "cell_type": "markdown", - "id": "d8b4cf19", - "metadata": {}, - "source": [ - "# NDRE Computation\n", - "In the same manner, we can compute, filter, and display another index useful for MicaSense cameras, 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, - "id": "b3fe7a83", - "metadata": {}, - "outputs": [], - "source": [ - "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n", - "rededge_band = thecapture.band_names_lower().index(\"red edge\")\n", - "ndre = (thelayer[:, :, nir_band] - thelayer[:, :, rededge_band]) / (\n", - " thelayer[:, :, nir_band] + thelayer[:, :, rededge_band]\n", - ")\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(np.percentile(masked_ndre, 0.5))\n", - "ndre_hist_max = np.max(np.percentile(masked_ndre, 99.5))\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(\n", - " gamma_corr_rgb,\n", - " masked_ndre,\n", - " figsize=(14, 7),\n", - " title=\"NDRE filtered to only plants over RGB base layer\",\n", - " vmin=min_display_ndre,\n", - " vmax=max_display_ndre,\n", - ")\n", - "fig.savefig(thecapture.uuid + \"_ndre_over_rgb.png\")" - ] - }, - { - "cell_type": "markdown", - "id": "898cc28a", - "metadata": {}, - "source": [ - "# Thermal Imagery\n", - "If our image is from an Altum or Altum-PT 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" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "176e2c34", - "metadata": {}, - "outputs": [], - "source": [ - "if len(thecapture.lw_indices()) > 0:\n", - " lwir_band = thecapture.band_names_lower().index(\"lwir\")\n", - "\n", - " # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n", - " if panchroCam:\n", - " masked_thermal = sharpened_stack[:, :, lwir_band]\n", - " else:\n", - " masked_thermal = im_aligned[:, :, lwir_band]\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", - " # Compute a histogram\n", - " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - " axis.hist(\n", - " masked_thermal.ravel(),\n", - " bins=512,\n", - " range=(\n", - " np.min(np.percentile(masked_thermal, 1)),\n", - " np.max(np.percentile(masked_thermal, 99)),\n", - " ),\n", - " )\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(\n", - " gamma_corr_rgb,\n", - " masked_thermal,\n", - " figsize=(14, 7),\n", - " title=\"Temperature over True Color\",\n", - " vmin=min_display_therm,\n", - " 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=0.4,\n", - " contour_fmt=\"%.0fC\",\n", - " )\n", - " fig.savefig(thecapture.uuid + \"_thermal_over_rgb.png\")" - ] - }, - { - "cell_type": "markdown", - "id": "bc9b29ec", - "metadata": {}, - "source": [ - "# Red vs NIR Reflectance\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, - "id": "c0250872", - "metadata": {}, - "outputs": [], - "source": [ - "x_band = red_band\n", - "y_band = nir_band\n", - "x_max = np.max(np.percentile(im_aligned[:, :, x_band], 99.99))\n", - "y_max = np.max(np.percentile(im_aligned[:, :, y_band], 99.99))\n", - "\n", - "fig = plt.figure(figsize=(12, 12))\n", - "plt.hexbin(\n", - " im_aligned[:, :, x_band],\n", - " im_aligned[:, :, y_band],\n", - " gridsize=640,\n", - " bins=\"log\",\n", - " extent=(0, x_max, 0, y_max),\n", - ")\n", - "ax = fig.gca()\n", - "ax.set_xlim([0, x_max])\n", - "ax.set_ylim([0, y_max])\n", - "plt.xlabel(\"{} Reflectance\".format(thecapture.band_names()[x_band]))\n", - "plt.ylabel(\"{} Reflectance\".format(thecapture.band_names()[y_band]))\n", - "plt.show()\n", - "\n", - "if panchroCam:\n", - " x_band = red_band\n", - " y_band = nir_band\n", - " x_max = np.max(np.percentile(sharpened_stack[:, :, x_band], 99.99))\n", - " y_max = np.max(np.percentile(sharpened_stack[:, :, y_band], 99.99))\n", - "\n", - " fig = plt.figure(figsize=(12, 12))\n", - " plt.hexbin(\n", - " sharpened_stack[:, :, x_band],\n", - " sharpened_stack[:, :, y_band],\n", - " gridsize=640,\n", - " bins=\"log\",\n", - " extent=(0, x_max, 0, y_max),\n", - " )\n", - " ax = fig.gca()\n", - " ax.set_xlim([0, x_max])\n", - " ax.set_ylim([0, y_max])\n", - " plt.xlabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[x_band]))\n", - " plt.ylabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[y_band]))\n", - " plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "4ec2d106", - "metadata": {}, - "source": [ - "# Print warp_matrices for usage elsewhere, such as Batch Processing\n", - "Lastly, 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, - "id": "c7cf8f48", - "metadata": {}, - "outputs": [], - "source": [ - "if panchroCam:\n", - " print(warp_matrices_SIFT)\n", - "else:\n", - " print(warp_matrices)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.12" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/Batch Processing v2.ipynb b/Batch Processing v2.ipynb index 7aa2ca1..969934f 100644 --- a/Batch Processing v2.ipynb +++ b/Batch Processing v2.ipynb @@ -1,384 +1,379 @@ { - "cells": [ - { - "cell_type": "markdown", - "id": "20e37d9b", - "metadata": {}, - "source": [ - "# Batch Processing Example\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. Part of this process (via `imageutils.write_exif_to_stack`) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.\n", - "\n", - "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support" - ] + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Batch Processing Example\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. Part of this process (via `imageutils.write_exif_to_stack`) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.\n", + "\n", + "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support" + ], + "id": "20e37d9b" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ], + "execution_count": null, + "outputs": [], + "id": "c64ead96" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load Images into ImageSet\n" + ], + "id": "80646655" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "from ipywidgets import FloatProgress, Layout\n", + "from IPython.display import display\n", + "import micasense.capture as capture\n", + "import os\n", + "from pathlib import Path\n", + "\n", + "# set to True if you have an Altum-PT\n", + "# or RedEdge-P and wish to output pan-sharpened stacks\n", + "panSharpen = True\n", + "\n", + "# If creating a lot of stacks, it is more efficient to save the metadata\n", + "# and then write all of the exif to the images after the stacks are created\n", + "write_exif_to_individual_stacks = False\n", + "\n", + "panelNames = None\n", + "useDLS = True\n", + "\n", + "# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html\n", + "imagePath = Path(\"./data/REDEDGE-MX\")\n", + "\n", + "# these will return lists of image paths as strings. Comment out of you aren't using panels.\n", + "panelNames = list(imagePath.glob(\"IMG_0001_*.tif\"))\n", + "panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "if panelNames:\n", + " panelCap = capture.Capture.from_filelist(panelNames)\n", + "\n", + "# destinations on your computer to put the stacks\n", + "# and RGB thumbnails\n", + "outputPath = imagePath / \"..\" / \"stacks\"\n", + "thumbnailPath = outputPath / \"thumbnails\"\n", + "\n", + "cam_model = panelCap.camera_model\n", + "cam_serial = panelCap.camera_serial\n", + "\n", + "# determine if this sensor has a panchromatic band\n", + "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", + " panchroCam = True\n", + "else:\n", + " panchroCam = False\n", + " panSharpen = False\n", + "\n", + "# if this is a multicamera system like the RedEdge-MX Dual,\n", + "# we can combine the two serial numbers to help identify\n", + "# this camera system later.\n", + "if len(panelCap.camera_serials) > 1:\n", + " cam_serial = \"_\".join(panelCap.camera_serials)\n", + " print(\"Serial number:\", cam_serial)\n", + "else:\n", + " cam_serial = panelCap.camera_serial\n", + " print(\"Serial number:\", cam_serial)\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; but 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(\n", + " v is None for v in panelCap.panel_albedo()\n", + " ):\n", + " panel_reflectance_by_band = panelCap.panel_albedo()\n", + " else:\n", + " panel_reflectance_by_band = [0.49] * len(\n", + " panelCap.eo_band_names()\n", + " ) # 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\"" + ], + "execution_count": null, + "outputs": [], + "id": "2081b06b" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import micasense.imageset as imageset\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", + "\n", + "\n", + "def update_f(val):\n", + " if (\n", + " val - f.value\n", + " ) > 0.005 or val == 1: # reduces cpu usage from updating the progressbar by 10x\n", + " f.value = val\n", + "\n", + "# time is not recognized and would remove the imageset as unused\n", + "imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", + "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", + "update_f(1.0)" + ], + "execution_count": null, + "outputs": [], + "id": "60736009" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Capture map\n", + "We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath." + ], + "id": "6c0a234d" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import numpy as np\n", + "from mapboxgl.viz import CircleViz\n", + "from mapboxgl.utils import df_to_geojson\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 access token here (https://account.mapbox.com/access-tokens/)\n", + "token = \"YOUR_MAPBOX_ACCESS_TOKEN\"\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", + "\n", + "geojson_data = df_to_geojson(df, columns[3:], lat=\"latitude\", lon=\"longitude\")\n", + "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", + "color_stops = create_color_stops(breaks, colors=\"YlOrRd\")\n", + "\n", + "viz = CircleViz(\n", + " geojson_data,\n", + " access_token=token,\n", + " color_property=color_property,\n", + " color_stops=color_stops,\n", + " center=[df[\"longitude\"].median(), df[\"latitude\"].median()],\n", + " zoom=16,\n", + " height=\"600px\",\n", + " style=\"mapbox://styles/mapbox/satellite-streets-v9\",\n", + ")\n", + "viz.show()" + ], + "execution_count": null, + "outputs": [], + "id": "9e9c437c" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Define which warp method to use\n", + "For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) 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. To use this simple alignment, simply set `warp_matrices=None` \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." + ], + "id": "e540a655" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "from micasense.warp_io import load_warp_matrices\n", + "\n", + "if panchroCam:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", + "else:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", + "\n", + "if Path(\"./\" + warp_matrices_filename).is_file():\n", + " print(\"Found existing warp matrices for camera\", cam_serial)\n", + " loaded = load_warp_matrices(\n", + " warp_matrices_filename, as_projective=panchroCam\n", + " )\n", + "\n", + " if panchroCam:\n", + " warp_matrices_SIFT = loaded\n", + " else:\n", + " warp_matrices = loaded\n", + " print(\"Loaded warp matrices from\", Path(\"./\" + warp_matrices_filename).resolve())\n", + "else:\n", + " print(\"No warp matrices found at expected location:\", warp_matrices_filename)" + ], + "execution_count": null, + "outputs": [], + "id": "635ad8d8" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Align images and save each capture to a layered TIFF file" + ], + "id": "3aae7f09" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import datetime\n", + "import micasense.imageutils as imageutils\n", + "\n", + "exif_list = []\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", + "\n", + "\n", + "def update_f2(val):\n", + " f2.value = val\n", + "\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 = str(i).zfill(4) + \"_\" + capture.uuid + \".tif\"\n", + " thumbnailFilename = str(i).zfill(4) + \"_\" + 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", + " if panchroCam:\n", + " capture.radiometric_pan_sharpened_aligned_capture(\n", + " warp_matrices=warp_matrices_SIFT,\n", + " irradiance_list=capture.dls_irradiace(),\n", + " img_type=img_type,\n", + " write_exif=write_exif_to_individual_stacks,\n", + " )\n", + " else:\n", + " capture.create_aligned_capture(\n", + " irradiance_list=irradiance, warp_matrices=warp_matrices\n", + " )\n", + " exif_list.append(\n", + " imageutils.prepare_exif_for_stacks(capture, fullOutputPath)\n", + " )\n", + " capture.save_capture_as_stack(\n", + " fullOutputPath,\n", + " pansharpen=panSharpen,\n", + " sort_by_wavelength=True,\n", + " write_exif=write_exif_to_individual_stacks,\n", + " )\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(\n", + " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", + " float(len(imgset.captures)) / float((end - start).total_seconds())\n", + " )\n", + ")" + ], + "execution_count": null, + "outputs": [], + "id": "cdc195da" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Write EXIF data to stacks\n", + "As mentioned above, it is more time intensive to write the exif data to each image as it is created. Here, we write the exif data after all of the TIFF files have been created. This should take a few seconds per stack." + ], + "id": "f4b88c80" + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if not write_exif_to_individual_stacks:\n", + " start = datetime.datetime.now()\n", + " for exif in exif_list:\n", + " imageutils.write_exif_to_stack(existing_exif_list=exif)\n", + " end = datetime.datetime.now()\n", + " print(\"Saving time: {}\".format(end - start))\n", + " print(\n", + " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", + " float(len(exif_list)) / float((end - start).total_seconds())\n", + " )\n", + " )" + ], + "execution_count": null, + "outputs": [], + "id": "fac88b6a" + } + ], + "metadata": { + "kernelspec": { + "display_name": "imageprocessing", + "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.12.12" + } }, - { - "cell_type": "code", - "execution_count": null, - "id": "c64ead96", - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "markdown", - "id": "80646655", - "metadata": {}, - "source": [ - "# Load Images into ImageSet\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2081b06b", - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import FloatProgress, Layout\n", - "from IPython.display import display\n", - "import micasense.capture as capture\n", - "import os\n", - "from pathlib import Path\n", - "\n", - "# set to True if you have an Altum-PT\n", - "# or RedEdge-P and wish to output pan-sharpened stacks\n", - "panSharpen = True\n", - "\n", - "# If creating a lot of stacks, it is more efficient to save the metadata\n", - "# and then write all of the exif to the images after the stacks are created\n", - "write_exif_to_individual_stacks = False\n", - "\n", - "panelNames = None\n", - "useDLS = True\n", - "\n", - "# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html\n", - "imagePath = Path(\"./data/REDEDGE-MX\")\n", - "\n", - "# these will return lists of image paths as strings. Comment out of you aren't using panels.\n", - "panelNames = list(imagePath.glob(\"IMG_0001_*.tif\"))\n", - "panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "if panelNames:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "\n", - "# destinations on your computer to put the stacks\n", - "# and RGB thumbnails\n", - "outputPath = imagePath / \"..\" / \"stacks\"\n", - "thumbnailPath = outputPath / \"thumbnails\"\n", - "\n", - "cam_model = panelCap.camera_model\n", - "cam_serial = panelCap.camera_serial\n", - "\n", - "# determine if this sensor has a panchromatic band\n", - "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", - " panchroCam = True\n", - "else:\n", - " panchroCam = False\n", - " panSharpen = False\n", - "\n", - "# if this is a multicamera system like the RedEdge-MX Dual,\n", - "# we can combine the two serial numbers to help identify\n", - "# this camera system later.\n", - "if len(panelCap.camera_serials) > 1:\n", - " cam_serial = \"_\".join(panelCap.camera_serials)\n", - " print(\"Serial number:\", cam_serial)\n", - "else:\n", - " cam_serial = panelCap.camera_serial\n", - " print(\"Serial number:\", cam_serial)\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; but 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(\n", - " v is None for v in panelCap.panel_albedo()\n", - " ):\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " panel_reflectance_by_band = [0.49] * len(\n", - " panelCap.eo_band_names()\n", - " ) # 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, - "id": "60736009", - "metadata": {}, - "outputs": [], - "source": [ - "import micasense.imageset as imageset\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", - "\n", - "\n", - "def update_f(val):\n", - " if (\n", - " val - f.value\n", - " ) > 0.005 or val == 1: # reduces cpu usage from updating the progressbar by 10x\n", - " f.value = val\n", - "\n", - "# time is not recognized and would remove the imageset as unused\n", - "imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "update_f(1.0)" - ] - }, - { - "cell_type": "markdown", - "id": "6c0a234d", - "metadata": {}, - "source": [ - "# Capture map\n", - "We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9e9c437c", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "from mapboxgl.viz import CircleViz\n", - "from mapboxgl.utils import df_to_geojson\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.eyJ1Ijoic3RlcGhlbm1hbmd1bTIiLCJhIjoiY2xmOXdnYzF1MDFqejNvdGE0YW13aTN5ZyJ9.AG_ckhUqTBjuGC2LuWCfQQ\"\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", - "\n", - "geojson_data = df_to_geojson(df, columns[3:], lat=\"latitude\", lon=\"longitude\")\n", - "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", - "color_stops = create_color_stops(breaks, colors=\"YlOrRd\")\n", - "\n", - "viz = CircleViz(\n", - " geojson_data,\n", - " access_token=token,\n", - " color_property=color_property,\n", - " color_stops=color_stops,\n", - " center=[df[\"longitude\"].median(), df[\"latitude\"].median()],\n", - " zoom=16,\n", - " height=\"600px\",\n", - " style=\"mapbox://styles/mapbox/satellite-streets-v9\",\n", - ")\n", - "viz.show()" - ] - }, - { - "cell_type": "markdown", - "id": "e540a655", - "metadata": {}, - "source": [ - "# Define which warp method to use\n", - "For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) 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. To use this simple alignment, simply set `warp_matrices=None` \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, - "id": "635ad8d8", - "metadata": {}, - "outputs": [], - "source": [ - "from skimage.transform import ProjectiveTransform\n", - "\n", - "if panchroCam:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", - "else:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", - "\n", - "if Path(\"./\" + warp_matrices_filename).is_file():\n", - " print(\"Found existing warp matrices for camera\", cam_serial)\n", - " load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True)\n", - " loaded_warp_matrices = []\n", - " for matrix in load_warp_matrices:\n", - " if panchroCam:\n", - " transform = ProjectiveTransform(matrix=matrix.astype(\"float64\"))\n", - " loaded_warp_matrices.append(transform)\n", - " else:\n", - " loaded_warp_matrices.append(matrix.astype(\"float32\"))\n", - "\n", - " if panchroCam:\n", - " warp_matrices_SIFT = loaded_warp_matrices\n", - " else:\n", - " warp_matrices = loaded_warp_matrices\n", - " print(\"Loaded warp matrices from\", Path(\"./\" + warp_matrices_filename).resolve())\n", - "else:\n", - " print(\"No warp matrices found at expected location:\", warp_matrices_filename)" - ] - }, - { - "cell_type": "markdown", - "id": "3aae7f09", - "metadata": {}, - "source": [ - "## Align images and save each capture to a layered TIFF file" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cdc195da", - "metadata": {}, - "outputs": [], - "source": [ - "import datetime\n", - "import micasense.imageutils as imageutils\n", - "\n", - "exif_list = []\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", - "\n", - "\n", - "def update_f2(val):\n", - " f2.value = val\n", - "\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 = str(i).zfill(4) + \"_\" + capture.uuid + \".tif\"\n", - " thumbnailFilename = str(i).zfill(4) + \"_\" + 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", - " if panchroCam:\n", - " capture.radiometric_pan_sharpened_aligned_capture(\n", - " warp_matrices=warp_matrices_SIFT,\n", - " irradiance_list=capture.dls_irradiace(),\n", - " img_type=img_type,\n", - " write_exif=write_exif_to_individual_stacks,\n", - " )\n", - " else:\n", - " capture.create_aligned_capture(\n", - " irradiance_list=irradiance, warp_matrices=warp_matrices\n", - " )\n", - " exif_list.append(\n", - " imageutils.prepare_exif_for_stacks(capture, fullOutputPath)\n", - " )\n", - " capture.save_capture_as_stack(\n", - " fullOutputPath,\n", - " pansharpen=panSharpen,\n", - " sort_by_wavelength=True,\n", - " write_exif=write_exif_to_individual_stacks,\n", - " )\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(\n", - " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", - " float(len(imgset.captures)) / float((end - start).total_seconds())\n", - " )\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "f4b88c80", - "metadata": {}, - "source": [ - "# Write EXIF data to stacks\n", - "As mentioned above, it is more time intensive to write the exif data to each image as it is created. Here, we write the exif data after all of the TIFF files have been created. This should take a few seconds per stack." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fac88b6a", - "metadata": {}, - "outputs": [], - "source": [ - "if not write_exif_to_individual_stacks:\n", - " start = datetime.datetime.now()\n", - " for exif in exif_list:\n", - " imageutils.write_exif_to_stack(existing_exif_list=exif)\n", - " end = datetime.datetime.now()\n", - " print(\"Saving time: {}\".format(end - start))\n", - " print(\n", - " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", - " float(len(exif_list)) / float((end - start).total_seconds())\n", - " )\n", - " )" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "imageprocessing", - "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.12.12" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 6c27f7f..d90aa78 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * fixed jupyter notebooks * `micasense.mp_config.spawn_pool()` using a per-pool spawn context (rawpy/OpenMP safe on Linux) * `micasense.warp_io` — `save_warp_matrices` / `load_warp_matrices` for `.npy` warp matrix I/O +* Alignment v2, Batch Processing v2 notebooks and `batch_processing_script.py` use `warp_io` ### Fixed diff --git a/batch_processing_script.py b/batch_processing_script.py index 0ea5c97..b43cd22 100644 --- a/batch_processing_script.py +++ b/batch_processing_script.py @@ -6,9 +6,8 @@ import numpy as np import pandas as pd from mapboxgl.utils import df_to_geojson -from skimage.transform import ProjectiveTransform - from micasense.capture import Capture +from micasense.warp_io import load_warp_matrices import micasense.imageset as imageset parser = argparse.ArgumentParser( @@ -114,19 +113,13 @@ if Path("./" + warp_matrices_filename).is_file(): print("Found existing warp matrices for camera", cam_serial) - load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True) - loaded_warp_matrices = [] - for matrix in load_warp_matrices: - if panchro_cam: - transform = ProjectiveTransform(matrix=matrix.astype("float64")) - loaded_warp_matrices.append(transform) - else: - loaded_warp_matrices.append(matrix.astype("float32")) - + loaded = load_warp_matrices( + warp_matrices_filename, as_projective=panchro_cam + ) if panchro_cam: - warp_matrices_SIFT = loaded_warp_matrices + warp_matrices_SIFT = loaded else: - warp_matrices = loaded_warp_matrices + warp_matrices = loaded print("Loaded warp matrices from", Path("./" + warp_matrices_filename).resolve()) else: print("No warp matrices found at expected location:", warp_matrices_filename) From 9d8f9d7506ff6eaa893799deee7dcb0fe3d8c158 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 15:15:18 +0000 Subject: [PATCH 4/9] Modernize test setup: pytest 9, coverage in CI, drop conda. Pin pytest>=9 with native pyproject config and coverage defaults, simplify CI to run pytest, and replace conda-based setup with venv/pip documentation. --- .github/workflows/publish.yml | 6 +- .github/workflows/python-test.yml | 6 +- CHANGELOG.md | 5 + MicaSense Image Processing Setup.ipynb | 134 ++++++++----------------- README.md | 14 +-- micasense_conda_env.yml | 28 ------ pyproject.toml | 20 +++- 7 files changed, 77 insertions(+), 136 deletions(-) delete mode 100644 micasense_conda_env.yml diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index 0a8c3ba..cf59f60 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -32,10 +32,10 @@ jobs: sudo apt-get -qq install exiftool libgdal-dev libzbar0t64 GDAL_VERSION=$(gdal-config --version) pip install "gdal==$GDAL_VERSION" + - name: Install package and test dependencies + run: pip install -e ".[test]" - name: Test with pytest - run: | - pip install -e .[test] - pytest -m "not newgdal" + run: pytest - name: Build package run: | python -m pip install build diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml index c3a0fa1..79a7c55 100644 --- a/.github/workflows/python-test.yml +++ b/.github/workflows/python-test.yml @@ -37,10 +37,10 @@ jobs: sudo apt-get -qq install exiftool libgdal-dev libzbar0t64 GDAL_VERSION=$(gdal-config --version) pip install "gdal==$GDAL_VERSION" + - name: Install package and test dependencies + run: pip install -e ".[test]" - name: Test with pytest - run: | - pip install -e .[test] - pytest -m "not newgdal" --cov=./ --cov-report=xml + run: pytest - uses: actions/upload-artifact@v6 if: ${{ matrix.python-version == '3.12' }} with: diff --git a/CHANGELOG.md b/CHANGELOG.md index d90aa78..5c21be6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * `micasense.mp_config.spawn_pool()` using a per-pool spawn context (rawpy/OpenMP safe on Linux) * `micasense.warp_io` — `save_warp_matrices` / `load_warp_matrices` for `.npy` warp matrix I/O * Alignment v2, Batch Processing v2 notebooks and `batch_processing_script.py` use `warp_io` +* pytest 9+ native `[tool.pytest]` config with coverage in CI + +### Changed + +* setup docs and README use venv/pip instead of conda; removed `micasense_conda_env.yml` ### Fixed diff --git a/MicaSense Image Processing Setup.ipynb b/MicaSense Image Processing Setup.ipynb index 74148d7..7341455 100644 --- a/MicaSense Image Processing Setup.ipynb +++ b/MicaSense Image Processing Setup.ipynb @@ -8,138 +8,86 @@ "\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", + "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.\n", + "Our tutorials are written using Python 3.10+. 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 its python wrapper to open and manipulate RedEdge images to transform raw digital number values into quantitative reflectance.\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 OSX, 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", + "The following software and libraries are required for this tutorial:\n", "\n", - "* [python3](https://www.python.org/downloads/release/latest)\n", + "* [python3](https://www.python.org/downloads/)\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/smarnach/pyexiftool) (version 0.4.13 only)\n", + "* [openCV](https://opencv.org/)\n", + "* [matplotlib](https://matplotlib.org/)\n", + "* [exiftool](https://exiftool.org/) + [pyexiftool](https://github.com/smarnach/pyexiftool)\n", "* [scikit-image](https://scikit-image.org/)\n", "* [zbar](https://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", + "* [GDAL](https://gdal.org/) (system package; match the pip `gdal` wheel to your installed version)\n", + "* [pandas](https://pandas.pydata.org/) (notebooks)\n", + "* [jupyter](https://jupyter.org/) (notebooks)\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", + "Install the library and notebook extras from a clone of this repository:\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", + " git clone https://github.com/brainergylab/micasense_imageprocessing.git\n", + " cd micasense_imageprocessing\n", + " python3 -m venv .venv\n", + " source .venv/bin/activate # Windows: .venv\\Scripts\\activate\n", + " pip install -e \".[test]\"\n", + " pip install jupyter pandas mapboxgl jenkspy imageio\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", + "On Linux, install system packages before `pip install` (see CI in `.github/workflows/python-test.yml`):\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", + " sudo apt-get update\n", + " sudo apt-get install exiftool libgdal-dev libzbar0\n", + " GDAL_VERSION=$(gdal-config --version)\n", + " pip install \"gdal==$GDAL_VERSION\"\n", "\n", - " wget https://exiftool.org/Image-ExifTool-12.57.tar.gz\n", - " tar -xvzf Image-ExifTool-12.57.tar.gz \n", - " cd Image-ExifTool-12.57/\n", - " perl Makefile.PL \n", - " make test\n", - " sudo make install\n", + "### Linux (Debian/Ubuntu)\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", + "The following is what we use on a fresh Ubuntu image. First install system tools and libraries:\n", "\n", - "Once these tools are installed, you can check out this repository and create the `micasense conda` environment:\n", + " sudo apt install git git-lfs libzbar0 libgdal-dev exiftool make\n", "\n", - " git clone https://github.com/micasense/imageprocessing.git\n", - " cd imageprocessing\n", - " conda env create -f micasense_conda_env.yml\n", + "Clone the repository, create a virtual environment, and install Python dependencies as shown above.\n", "\n", - "Finally, one way to verify our install by running the built in tests:\n", + "Verify the install:\n", "\n", - " cd imageprocessing\n", - " conda activate micasense\n", - " pytest .\n", + " pytest\n", "\n", - "Or, to start working with the notebooks (including running the test code below):\n", + "Or start the notebooks:\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.7** version\n", + "* Install [Python 3](https://www.python.org/downloads/) for your system.\n", + "* Download the [exiftool windows package](https://exiftool.org/) and unzip it to a permanent location such as `c:\\exiftool\\`.\n", + "* Add an environment variable `exiftoolpath` pointing to `c:\\exiftool\\exiftool.exe` (Settings \u2192 Environment Variables \u2192 New).\n", + "* Open Command Prompt or PowerShell, `cd` to the repository, create and activate a venv, then run `pip install -e \".[test]\"` and install GDAL matching your environment if needed.\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", - "The following steps to get going on MacOS worked for us. \n", - "\n", - "First we installed `git` by installing the `xcode` developer tools, or you can follow the instructions at the [git site](https://git-scm.com/downloads). \n", - "\n", - "Next, we [downloaded and installed exiftool](https://exiftool.org/) using the MacOS installer.\n", - "\n", - "Third, we installed [Homebrew](https://brew.sh/) and used it to install `zbar`:\n", - "\n", - " /usr/bin/ruby -e \"$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)\"\n", - " brew install zbar\n", - " \n", - "Then we installed miniconda. If you're comfortable on the command line, navigate to the [miniconda download page](https://conda.io/miniconda.html) and download the installer for your system. Open an iTerm and follow the [installation instructions](https://conda.io/docs/user-guide/install/index.html).\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 an iTerm 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", + "* Install [git](https://git-scm.com/downloads) and [git-lfs](https://git-lfs.github.com/).\n", + "* Install [exiftool](https://exiftool.org/).\n", + "* Install [Homebrew](https://brew.sh/) and `zbar` / GDAL as needed: `brew install zbar gdal`\n", + "* Create a virtual environment in the repository and `pip install -e \".[test]\"` as above.\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", + "* Activate your virtual environment (`source .venv/bin/activate` or `.venv\\Scripts\\activate` on Windows)\n", + "* `cd` to the repository 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" + "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", + "\n" ] }, { diff --git a/README.md b/README.md index 1e92da7..5974adc 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ This package is available on PyPi and can be installed using `pip install micase The code in these tutorials consists of two parts. First, the tutorials generally end in `.ipynb` and are the Jupyter notebooks that were used to create the web page tutorials linked above. You can run this code by opening a -terminal/iTerm (linux/mac) or Anaconda Command Prompt (Windows), navigating to the folder you cloned the git repository +terminal (linux/mac) or Command Prompt/PowerShell (Windows), navigating to the folder you cloned the git repository into, and running ```bash @@ -95,8 +95,8 @@ Want to correct an issue or expand library functionality? Fork the repository, m on github. Have a question? Please double-check that you're able to run the setup notebook successfully, and resolve any issues -with that first. If you're pulling newer code, it may be necessary in some cases to delete and re-create -your `micasense` conda environment to make sure you have all of the expected packages. +with that first. If you're pulling newer code, recreate your virtual environment or reinstall with +`pip install -e ".[test]"` to pick up dependency changes. This code is a community effort and is not supported by MicaSense support. Please don't reach out to MicaSense support for issues with this codebase; instead, work through the above troubleshooting steps and @@ -104,15 +104,17 @@ then [create an issue on github](https://github.com/brainergylab/micasense_image ### Tests -Tests for many library functions are included in the `tests` diretory. Install the `pytest` module through your package -manager (e.g. `pip install pytest`) and then tests can be run from the main directory using the command: +Tests for many library functions are included in the `tests` directory. Install test dependencies and run from the repository root: ```bash +pip install -e ".[test]" pytest ``` +Coverage and default options are configured in `pyproject.toml` (`pytest>=9`). + Test execution can be relatively slow (2-3 minutes) as there is a lot of image processing occuring in some of the tests, -and quite a bit of re-used IO. To speed up tests, install the `pytest-xdist` plugin using `conda` or `pip` and achieve a +and quite a bit of re-used IO. To speed up tests, install the `pytest-xdist` plugin using `pip` and achieve a significant speed up by running tests in parallel. ```bash diff --git a/micasense_conda_env.yml b/micasense_conda_env.yml deleted file mode 100644 index 845e1f7..0000000 --- a/micasense_conda_env.yml +++ /dev/null @@ -1,28 +0,0 @@ -name: micasense -channels: - - conda-forge -dependencies: - - python=3.8 - - opencv=4.8 - - numpy - - jupyter - - matplotlib - - pysolar - - gdal - - scikit-image - - pytest - - git - - git-lfs - - pandas - - imageio - - nb_conda - - requests - - packaging - - pip - - pytest-xdist - - pip: - - pyzbar - - mapboxgl - - pyexiftool - - jenkspy - - rawpy diff --git a/pyproject.toml b/pyproject.toml index b506f16..2ae6db3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,8 +47,8 @@ Homepage = "https://github.com/brainergylab/micasense_imageprocessing" [project.optional-dependencies] test = [ - "pytest >= 7.4.3", - "pytest-cov", + "pytest>=9", + "pytest-cov>=6", ] docs = [ "jupyter", @@ -65,5 +65,19 @@ exclude = ["data*", "tests*"] [tool.pytest] testpaths = ["tests"] markers = [ - "newgdal", + "newgdal: tests that require newer GDAL stack behavior", ] +addopts = [ + "-m", + "not newgdal", + "--cov=micasense", + "--cov-report=xml", + "--cov-report=term-missing", +] + +[tool.coverage.run] +source = ["micasense"] + +[tool.coverage.report] +show_missing = true +skip_covered = true From c3396307b3f880b58bcbb29015a5a4b48a21ad84 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 19:07:08 +0000 Subject: [PATCH 5/9] Remove redundant work in SIFT_align_capture. Drop unused lists, duplicate raw() calls, and per-band undistort reloads without changing alignment results. --- CHANGELOG.md | 1 + micasense/capture.py | 25 ++++++++++--------------- 2 files changed, 11 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c21be6..947a97a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * `ImageSet.save_stacks` and `imageutils` alignment pools use spawn context (not global `set_start_method`) * declare `rawpy` as a runtime dependency * remove unused `Capture.__sift_warp_matrices` +* `SIFT_align_capture`: remove redundant raw/undistort work and dead accumulators ## [0.1.1] - 2025-12-28 diff --git a/micasense/capture.py b/micasense/capture.py index 0f5bdb6..88e85e7 100644 --- a/micasense/capture.py +++ b/micasense/capture.py @@ -1159,14 +1159,16 @@ def SIFT_align_capture( descriptors = [] img_index = list(range(len(self.images))) img_index.pop(ref) - ref_shape = self.images[ref].raw().shape + ref_image = self.images[ref] + ref_raw = ref_image.raw() + ref_shape = ref_raw.shape rest_shape = self.images[img_index[0]].raw().shape scale = np.array(ref_shape) / np.array(rest_shape) # use the calibrated warp matrices to verify keypoints warp_matrices_calibrated = self.get_warp_matrices(ref_index=ref) - ref_img = self.images[ref].undistorted(self.images[ref].raw()) + ref_img = ref_image.undistorted(ref_raw) if ref_shape != rest_shape: ref_img = resize(ref_img, rest_shape) ref_image_SIFT = (ref_img / ref_img.max() * 65535).astype(np.uint16) @@ -1176,16 +1178,18 @@ def SIFT_align_capture( descriptor_ref = descriptor_extractor.descriptors if verbose > 1: logger.info("found %d keypoints in the reference image", len(keypoints_ref)) - match_images = [] ratio = [] filter_tr = [] + raw_shapes = {} img_index = np.array(img_index) # extract keypoints % descriptors for ix in img_index: - img = self.images[ix].undistorted(self.images[ix].raw()) + raw_band = self.images[ix].raw() + raw_shapes[ix] = raw_band.shape + img = self.images[ix].undistorted(raw_band) if not img.shape == rest_shape: # if we have a thermal image, upsample to match the resolution of the multispec images - img_base = self.images[ix].raw()[self.images[ix].raw() > 0].min() + img_base = raw_band[raw_band > 0].min() img = img.astype(float) img[img > 0] = img[img > 0] - img_base img = resize(img, rest_shape) @@ -1199,7 +1203,6 @@ def SIFT_align_capture( else: # less strict filtering for the BLUE images filter_tr.append(err_blue) - match_images.append(img) descriptor_extractor.detect_and_extract(img) keypoints.append(descriptor_extractor.keypoints) descriptors.append(descriptor_extractor.descriptors) @@ -1248,11 +1251,9 @@ def SIFT_align_capture( warp_matrices_calibrated[ix] = np.dot(warp_blue_ref, warpBLUE[ix]) models = [] - kp_image = [] - kp_ref = [] for m, k, ix, t in zip(matches, keypoints, img_index, filter_tr): # we need to down scale the thermal image for the proper transform - scale_i = np.array(self.images[ix].raw().shape) / np.array(rest_shape) + scale_i = np.array(raw_shapes[ix]) / np.array(rest_shape) filtered_kpi, filtered_kpr, filtered_match, err = self.filter_keypoints( k, @@ -1286,13 +1287,7 @@ def SIFT_align_capture( ix, self.images[ix].band_name, ) - kpi, kpr = filtered_kpi, filtered_kpr models.append(P) - kp_image.append(kpi) - kp_ref.append(kpr) - img = self.images[ix].undistorted(self.images[ix].raw()) - - # no need for the upsampled stacks here if verbose > 0: logger.info("Finished aligning band %d", ix) From cbca97926574b4b883788693a0b65ac1447e6d11 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Wed, 3 Jun 2026 19:07:47 +0000 Subject: [PATCH 6/9] Ignore local test and venv artifacts in git. Add .coverage, coverage.xml, and .venv/ to .gitignore. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 6843774..7a6a429 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,9 @@ # Add any directories, files, or patterns you don't want to be tracked by version control *.pyc .ipynb_checkpoints +.coverage +coverage.xml +.venv/ __pycache__ /dist/ /*.egg-info From 91b0f8715f032efca4d2e6e90d26eca39dd4f0d8 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Fri, 19 Jun 2026 08:09:47 +0200 Subject: [PATCH 7/9] update and run pre-commit --- .pre-commit-config.yaml | 6 +- Alignment v2.ipynb | 1876 ++++++++++++++++++------------------ Batch Processing v2.ipynb | 753 ++++++++------- batch_processing_script.py | 5 +- micasense/imageutils.py | 1 - micasense/warp_io.py | 8 +- tests/test_sift_align.py | 2 - tests/test_warp_io.py | 1 - 8 files changed, 1319 insertions(+), 1333 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index efb4d41..27c72c2 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,13 +1,13 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v5.0.0 + rev: v6.0.0 hooks: - id: end-of-file-fixer - id: trailing-whitespace - id: check-illegal-windows-names - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.14.6 + rev: v0.15.18 hooks: # Run the linter. - id: ruff @@ -17,7 +17,7 @@ repos: - id: ruff-format types_or: [ python, pyi ] # avoid jupyter - repo: https://github.com/codespell-project/codespell - rev: v2.3.0 + rev: v2.4.2 hooks: - id: codespell additional_dependencies: diff --git a/Alignment v2.ipynb b/Alignment v2.ipynb index 88676ac..93d18d9 100644 --- a/Alignment v2.ipynb +++ b/Alignment v2.ipynb @@ -1,941 +1,939 @@ { - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Active Image Alignment\n", - "For 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 RedEdge (3,M,MX) versions above 3.2.0 which include the \"RigRelatives\" tags, and all RedEdge-P/Altum/Altum-PT imagery. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.\n", - "\n", - "# Opening Images\n", - "As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize all the bands of a MicaSense 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 the `autoreload` extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time." - ], - "id": "573422ea" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ], - "execution_count": null, - "outputs": [], - "id": "b09cd3b1" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import micasense.capture as capture\n", - "\n", - "%matplotlib inline\n", - "from pathlib import Path\n", - "import matplotlib.pyplot as plt\n", - "\n", - "plt.rcParams[\"figure.facecolor\"] = \"w\"\n", - "\n", - "panelNames = None\n", - "\n", - "# set your image paths here. See more here: https://docs.python.org/3/library/pathlib.html\n", - "# if using Windows, you need to an an \"r\" to the path like this: Path(r\"C:\\Files\")\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-MX-DUAL\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0431_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-MX\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0020_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0001_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/ALTUM\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0021_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "# imagePath = Path(\"./data/REDEDGE-P\")\n", - "\n", - "# # these will return lists of image paths as strings\n", - "# imageNames = list(imagePath.glob('IMG_0011_*.tif'))\n", - "# imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", - "# panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "imagePath = Path(\"./data/ALTUM-PT\")\n", - "\n", - "# these will return lists of image paths as strings\n", - "imageNames = list(imagePath.glob(\"IMG_0010_*.tif\"))\n", - "imageNames = [x.as_posix() for x in imageNames]\n", - "\n", - "panelNames = list(imagePath.glob(\"IMG_0000_*.tif\"))\n", - "panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "\n", - "if panelNames is not None:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "else:\n", - " panelCap = None\n", - "\n", - "thecapture = capture.Capture.from_filelist(imageNames)\n", - "\n", - "# get camera model for future use\n", - "cam_model = thecapture.camera_model\n", - "# if this is a multicamera system like the RedEdge-MX Dual,\n", - "# we can combine the two serial numbers to help identify\n", - "# this camera system later.\n", - "if len(thecapture.camera_serials) > 1:\n", - " cam_serial = \"_\".join(thecapture.camera_serials)\n", - " print(cam_serial)\n", - "else:\n", - " cam_serial = thecapture.camera_serial\n", - "\n", - "print(\"Camera model:\", cam_model)\n", - "print(\"Bit depth:\", thecapture.bits_per_pixel)\n", - "print(\"Camera serial number:\", cam_serial)\n", - "print(\"Capture ID:\", thecapture.uuid)\n", - "\n", - "# determine if this sensor has a panchromatic band\n", - "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", - " panchroCam = True\n", - "else:\n", - " panchroCam = False\n", - " panSharpen = False\n", - "\n", - "if panelCap is not None:\n", - " if panelCap.panel_albedo() is not None:\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " panel_reflectance_by_band = [0.49] * len(\n", - " thecapture.eo_band_names()\n", - " ) # RedEdge band_index order\n", - " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)\n", - " irradiance_list = panelCap.panel_irradiance(panel_reflectance_by_band) + [\n", - " 0\n", - " ] # add to account for uncalibrated LWIR band, if applicable\n", - " img_type = \"reflectance\"\n", - " thecapture.plot_undistorted_reflectance(panel_irradiance)\n", - "else:\n", - " if thecapture.dls_present():\n", - " img_type = \"reflectance\"\n", - " irradiance_list = thecapture.dls_irradiance() + [0]\n", - " thecapture.plot_undistorted_reflectance(thecapture.dls_irradiance())\n", - " else:\n", - " img_type = \"radiance\"\n", - " thecapture.plot_undistorted_radiance()\n", - " irradiance_list = None" - ], - "execution_count": null, - "outputs": [], - "id": "7426068b" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Check for existing warp matrices \n", - "If we have already successfully aligned captures from this specific camera, we can typically save some time and use the alignment warp matrices for other captures from the same camera" - ], - "id": "8de5a774" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "from micasense.warp_io import load_warp_matrices\n", - "\n", - "if panchroCam:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", - "else:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", - "\n", - "if Path(\"./\" + warp_matrices_filename).is_file():\n", - " print(\"Found existing warp matrices for camera\", cam_serial)\n", - " loaded = load_warp_matrices(\n", - " warp_matrices_filename, as_projective=panchroCam\n", - " )\n", - " print(\"Warp matrices successfully loaded.\")\n", - "\n", - " if panchroCam:\n", - " warp_matrices_SIFT = loaded\n", - " else:\n", - " warp_matrices = loaded\n", - "else:\n", - " print(\"No existing warp matrices found. Create them later in the notebook.\")\n", - " warp_matrices_SIFT = False\n", - " warp_matrices = False" - ], - "execution_count": null, - "outputs": [], - "id": "35de1be0" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Unwarp and Align (OpenCV method for RedEdge3/M/MX/Dual and original Altum)\n", - "Alignment is a three-step process:\n", - "\n", - "Images are unwarped using the built-in lens calibration\n", - "A transformation is found to align each band to a common band\n", - "The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", - "We provide utilities to find the alignment transformations within a single capture. Our experience shows that once a good alignment 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 alignment 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." - ], - "id": "185627ff" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import cv2\n", - "import time\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", - "# We use a different alignment method for RedEdge-P and Altum-PT,\n", - "# so if the imagery is from this kind of camera, skip this step\n", - "if not panchroCam:\n", - " st = time.time()\n", - " # set to True if you'd like to ignore existing warp matrices and create new ones\n", - " regenerate = True\n", - " pyramid_levels = (\n", - " 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment\n", - " )\n", - " max_alignment_iterations = 10\n", - "\n", - " # match_index:\n", - " # for non-panchromatic cameras we want to use band 1, which is green\n", - " # the green band has zero rig relative offsets\n", - " # NOTE: These band numbers are zero-indexed\n", - " # special parameters for RedEdge-MX dual camera system\n", - " if len(thecapture.eo_band_names()) == 10:\n", - " print(\"is 10 band\")\n", - " pyramid_levels = 3\n", - " match_index = 4\n", - " max_alignment_iterations = 20\n", - " else:\n", - " match_index = 1\n", - "\n", - " warp_mode = (\n", - " cv2.MOTION_HOMOGRAPHY\n", - " ) # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", - "\n", - " print(\n", - " \"Aligning images. Depending on settings this can take from a few seconds to many minutes\"\n", - " )\n", - " # Can potentially increase max_iterations for better results, but longer runtimes\n", - " if warp_matrices and not regenerate:\n", - " print(\"Using existing warp matrices...\")\n", - " try:\n", - " irradiance = panel_irradiance + [0]\n", - " except NameError:\n", - " irradiance = None\n", - " else:\n", - " warp_matrices, alignment_pairs = imageutils.align_capture(\n", - " thecapture,\n", - " ref_index=match_index,\n", - " max_iterations=max_alignment_iterations,\n", - " warp_mode=warp_mode,\n", - " pyramid_levels=pyramid_levels,\n", - " )\n", - "\n", - " print(\"Finished Aligning\")\n", - " et = time.time()\n", - " elapsed_time = et - st\n", - " print(\"Alignment time:\", int(elapsed_time), \"seconds\")" - ], - "execution_count": null, - "outputs": [], - "id": "e154cdd1" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Crop Aligned Images and create aligned capture stack \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 take 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 to 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." - ], - "id": "c06bfdf3" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "if not panchroCam:\n", - " cropped_dimensions, edges = imageutils.find_crop_bounds(\n", - " thecapture, warp_matrices, warp_mode=warp_mode, reference_band=match_index\n", - " )\n", - " print(\"Cropped dimensions:\", cropped_dimensions)\n", - " im_aligned = thecapture.create_aligned_capture(\n", - " warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type\n", - " )" - ], - "execution_count": null, - "outputs": [], - "id": "5fbfb346" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Alignment and pan-sharpening for RedEdge-P and Altum-PT\n", - "For older cameras we use OpenCV for capture alignment. For RedEdge-P and Altum-PT, we use SIFT (scale-invariant feature transform). We will use SIFT to create the warp matrices, then use these warp matrices to align the capture. See more here: \n", - "https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_sift.html\n", - "https://en.wikipedia.org/wiki/Scale-invariant_feature_transform\n", - "\n", - "For sensors with a panchromatic band (Altum-PT or RedEdge-P), we may wish to create a pan-sharpened stack. This example uses a linear interpolation method for pan-sharpening.\n", - "\n", - "The `radiometric_pan_sharpen` function takes in the SIFT warp matrices and outputs an upsampled image stack (all bands are changed to the resolution of the panchromatic sensor) and a pan-sharpened stack. \n", - "\n", - "Note: it is important to choose an initial alignment image that has a lot of straight, manmade features, such as roads or buildings. Trying to align a capture that is strictly vegetation will take a long time and may not be very good. This process can take a while, depending on how feature-rich the capture is. We have seen some Altum-PT captures align after 2 minutes, and some can take upwards of 30 minutes. Be patient! Once an alignment has been found, it can be used on other captures from the same camera, and the alignment process will be much faster. " - ], - "id": "c362623a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# from micasense.imageutils import brovey_pan_sharpen,radiometric_pan_sharpen\n", - "from skimage.transform import ProjectiveTransform\n", - "import time\n", - "\n", - "if panchroCam:\n", - " # set to True if you'd like to ignore existing warp matrices and create new ones\n", - " regenerate = True\n", - " st = time.time()\n", - " if not warp_matrices_SIFT or regenerate:\n", - " print(\"Generating new warp matrices...\")\n", - " warp_matrices_SIFT = thecapture.SIFT_align_capture(min_matches=10)\n", - "\n", - " sharpened_stack, upsampled = thecapture.radiometric_pan_sharpened_aligned_capture(\n", - " warp_matrices=warp_matrices_SIFT,\n", - " irradiance_list=irradiance_list,\n", - " img_type=img_type,\n", - " )\n", - "\n", - " # we can also use the Rig Relatives from the image metadata to do a quick, rudimentary alignment\n", - " # warp_matrices0=thecapture.get_warp_matrices(ref_index=5)\n", - " # sharpened_stack,upsampled = radiometric_pan_sharpen(thecapture,warp_matrices=warp_matrices0)\n", - "\n", - " print(\"Pansharpened shape:\", sharpened_stack.shape)\n", - " print(\"Upsampled shape:\", upsampled.shape)\n", - " # re-assign to im_aligned to match rest of code\n", - " im_aligned = upsampled\n", - " et = time.time()\n", - " elapsed_time = et - st\n", - " print(\"Alignment and pan-sharpening time:\", int(elapsed_time), \"seconds\")" - ], - "execution_count": null, - "outputs": [], - "id": "12c2a278" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Save warp matrices\n", - "Once an alignment for your camera has been found, it can be saved to a file for later use with this notebook. It can also be used on the Batch Alignment notebook for aligning all of the captures from an entire flight." - ], - "id": "3e9d51b1" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "from micasense.warp_io import save_warp_matrices\n", - "\n", - "if panchroCam:\n", - " working_wm = warp_matrices_SIFT\n", - "else:\n", - " working_wm = warp_matrices\n", - "if not Path(\"./\" + warp_matrices_filename).is_file() or regenerate:\n", - " save_warp_matrices(warp_matrices_filename, working_wm)\n", - " print(\"Saved to\", Path(\"./\" + warp_matrices_filename).resolve())\n", - "else:\n", - " print(\"Matrices already exist at\", Path(\"./\" + warp_matrices_filename).resolve())" - ], - "execution_count": null, - "outputs": [], - "id": "84855db4" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Multispectral band histogram \n", - "We can compare the radiance between multispectral bands, and between upsampled/pansharpened in the case of RedEdge-P and Altum-PT" - ], - "id": "47d3fb77" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "theColors = {\n", - " \"Blue\": \"blue\",\n", - " \"Green\": \"green\",\n", - " \"Red\": \"red\",\n", - " \"Red edge\": \"maroon\",\n", - " \"NIR\": \"purple\",\n", - " \"Panchro\": \"yellow\",\n", - " \"PanchroB\": \"orange\",\n", - " \"Red edge-740\": \"salmon\",\n", - " \"Red Edge\": \"maroon\",\n", - " \"Blue-444\": \"aqua\",\n", - " \"Green-531\": \"lime\",\n", - " \"Red-650\": \"lightcoral\",\n", - " \"Red edge-705\": \"brown\",\n", - "}\n", - "\n", - "eo_count = len(thecapture.eo_indices())\n", - "multispec_min = np.min(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 0.01))\n", - "multispec_max = np.max(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 99.99))\n", - "\n", - "theRange = (multispec_min, multispec_max)\n", - "\n", - "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - "for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", - " axis.hist(\n", - " im_aligned[:, :, x].ravel(),\n", - " bins=512,\n", - " range=theRange,\n", - " histtype=\"step\",\n", - " label=y,\n", - " color=theColors[y],\n", - " linewidth=1.5,\n", - " )\n", - "plt.title(\"Multispectral histogram (radiance)\")\n", - "axis.legend()\n", - "plt.show()\n", - "\n", - "if panchroCam:\n", - " eo_count = len(thecapture.eo_indices())\n", - " multispec_min = np.min(\n", - " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 0.01)\n", - " )\n", - " multispec_max = np.max(\n", - " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 99.99)\n", - " )\n", - "\n", - " theRange = (multispec_min, multispec_max)\n", - "\n", - " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - " for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", - " axis.hist(\n", - " sharpened_stack[:, :, x].ravel(),\n", - " bins=512,\n", - " range=theRange,\n", - " histtype=\"step\",\n", - " label=y,\n", - " color=theColors[y],\n", - " linewidth=1.5,\n", - " )\n", - " plt.title(\"Pan-sharpened multispectral histogram (radiance)\")\n", - " axis.legend()\n", - " plt.show()" - ], - "execution_count": null, - "outputs": [], - "id": "8e036b45" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Visualize Aligned Images\n", - "Once the transformation has been found, it can be verified by compositing 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 external software such as QGIS. Useful composites are a naturally colored RGB as well as color infrared, or CIR." - ], - "id": "3ccb819b" - }, - { - "cell_type": "code", - "metadata": {}, - "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 = [\n", - " thecapture.band_names_lower().index(\"red\"),\n", - " thecapture.band_names_lower().index(\"green\"),\n", - " thecapture.band_names_lower().index(\"blue\"),\n", - "]\n", - "cir_band_indices = [\n", - " thecapture.band_names_lower().index(\"nir\"),\n", - " thecapture.band_names_lower().index(\"red\"),\n", - " thecapture.band_names_lower().index(\"green\"),\n", - "]\n", - "\n", - "# Create normalized stacks for viewing\n", - "im_display = np.zeros(\n", - " (im_aligned.shape[0], im_aligned.shape[1], im_aligned.shape[2]), dtype=float\n", - ")\n", - "im_min = np.percentile(\n", - " im_aligned[:, :, rgb_band_indices].flatten(), 0.5\n", - ") # modify these percentiles to adjust contrast\n", - "im_max = np.percentile(\n", - " im_aligned[:, :, rgb_band_indices].flatten(), 99.5\n", - ") # for many images, 0.5 and 99.5 are good values\n", - "\n", - "if panchroCam:\n", - " im_display_sharp = np.zeros(\n", - " (sharpened_stack.shape[0], sharpened_stack.shape[1], sharpened_stack.shape[2]),\n", - " dtype=float,\n", - " )\n", - " im_min_sharp = np.percentile(\n", - " sharpened_stack[:, :, rgb_band_indices].flatten(), 0.5\n", - " ) # modify these percentiles to adjust contrast\n", - " im_max_sharp = np.percentile(\n", - " sharpened_stack[:, :, rgb_band_indices].flatten(), 99.5\n", - " ) # for many images, 0.5 and 99.5 are good values\n", - "\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", - " if panchroCam:\n", - " im_display_sharp[:, :, i] = imageutils.normalize(\n", - " sharpened_stack[:, :, i], im_min_sharp, im_max_sharp\n", - " )\n", - "\n", - "rgb = im_display[:, :, rgb_band_indices]\n", - "\n", - "if panchroCam:\n", - " rgb_sharp = im_display_sharp[:, :, rgb_band_indices]\n", - "\n", - "nir_band = thecapture.band_names_lower().index(\"nir\")\n", - "red_band = thecapture.band_names_lower().index(\"red\")\n", - "\n", - "ndvi = (im_aligned[:, :, nir_band] - im_aligned[:, :, red_band]) / (\n", - " im_aligned[:, :, nir_band] + im_aligned[:, :, red_band]\n", - ")\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", - "if panchroCam:\n", - " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)\n", - "else:\n", - " fig, ax1 = plt.subplots(1, 1, figsize=figsize)\n", - "ax1.set_title(\"Red-Green-Blue Composite\")\n", - "ax1.imshow(rgb)\n", - "if panchroCam:\n", - " ax2.set_title(\"Red-Green-Blue Composite (pan-sharpened)\")\n", - " ax2.imshow(rgb_sharp)\n", - "\n", - "fig, (ax3, ax4) = plt.subplots(1, 2, figsize=figsize)\n", - "ax3.set_title(\"NDVI\")\n", - "ax3.imshow(ndvi)\n", - "ax4.set_title(\"Color Infrared (CIR) Composite\")\n", - "ax4.imshow(cir)\n", - "\n", - "# set custom lims if you want to zoom in to image to see more detail\n", - "# this is useful for comparing upsampled and pan-sharpened stacks\n", - "# custom_xlim=(1500,2000)\n", - "# custom_ylim=(2000,1500)\n", - "# plt.setp([ax1,ax2,ax3,ax4], xlim=custom_xlim, ylim=custom_ylim)\n", - "\n", - "plt.show()" - ], - "execution_count": null, - "outputs": [], - "id": "8626c90b" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Image Enhancement\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." - ], - "id": "17f71519" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "if panchroCam:\n", - " rgb = rgb_sharp\n", - "# 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()" - ], - "execution_count": null, - "outputs": [], - "id": "1185fd16" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Stack Export\n", - "We can easily export the stacks into an image 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. The stacks include geographic information. \n", - "\n", - "If you prefer, you may set `sort_by_wavelength` to `True` in the `save_capture_as_stack` function.\n", - "\n", - "Unless otherwise specified, this will save in your working folder, that is your `imageprocessing` directory." - ], - "id": "6d1f6189" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# set output name to unique capture ID, e.g. FWoNSvgDNBX63Xv378qs\n", - "outputName = thecapture.uuid\n", - "\n", - "st = time.time()\n", - "if panchroCam:\n", - " # in this example, we can export both a pan-sharpened stack and an upsampled stack\n", - " # so you can compare them in GIS. In practice, you would typically only output the pansharpened stack\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-pansharpened.tif\", sort_by_wavelength=True, pansharpen=True\n", - " )\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-upsampled.tif\", sort_by_wavelength=True, pansharpen=False\n", - " )\n", - "else:\n", - " thecapture.save_capture_as_stack(\n", - " outputName + \"-noPanels.tif\", sort_by_wavelength=True\n", - " )\n", - "\n", - "et = time.time()\n", - "elapsed_time = et - st\n", - "print(\"Time to save stacks:\", int(elapsed_time), \"seconds.\")" - ], - "execution_count": null, - "outputs": [], - "id": "a2d716be" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# NDVI Computation\n", - "For raw index computation on single images, the `numpy` package provides a simple way to do math and simple visualization on images. Below, we compute and visualize an image histogram, and then use that to pick a color map range for visualizing the NDVI of an image.\n", - "\n", - "## Plant Classification\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." - ], - "id": "f4d9cd64" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "nir_band = thecapture.band_names_lower().index(\"nir\")\n", - "red_band = thecapture.band_names_lower().index(\"red\")\n", - "\n", - "thelayer = im_aligned\n", - "if panchroCam:\n", - " thelayer = sharpened_stack\n", - "np.seterr(\n", - " divide=\"ignore\", invalid=\"ignore\"\n", - ") # 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 = (thelayer[:, :, nir_band] - thelayer[:, :, red_band]) / (\n", - " thelayer[:, :, nir_band] + thelayer[:, :, red_band]\n", - ")\n", - "print(\"Image type:\", img_type)\n", - "\n", - "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", - "# this does not seem to work on panchro stacks\n", - "if img_type == \"reflectance\":\n", - " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < 0.20, ndvi)\n", - "elif img_type == \"radiance\":\n", - " lower_pct_radiance = np.percentile(thelayer[:, :, nir_band], 10.0)\n", - " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < lower_pct_radiance, ndvi)\n", - "# Compute and display a histogram\n", - "# ndvi_hist_min = np.min(ndvi)\n", - "# ndvi_hist_max = np.max(ndvi)\n", - "ndvi_hist_min = np.min(np.percentile(ndvi, 0.5))\n", - "ndvi_hist_max = np.max(np.percentile(ndvi, 99.5))\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(\n", - " ndvi.flatten(), 99.5\n", - ") # 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(\n", - " gamma_corr_rgb,\n", - " masked_ndvi,\n", - " figsize=(14, 7),\n", - " title=\"NDVI filtered to only plants over RGB base layer\",\n", - " vmin=min_display_ndvi,\n", - " vmax=max_display_ndvi,\n", - ")\n", - "fig.savefig(thecapture.uuid + \"_ndvi_over_rgb.png\")" - ], - "execution_count": null, - "outputs": [], - "id": "4a572a55" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# NDRE Computation\n", - "In the same manner, we can compute, filter, and display another index useful for MicaSense cameras, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health." - ], - "id": "d8b4cf19" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n", - "rededge_band = thecapture.band_names_lower().index(\"red edge\")\n", - "ndre = (thelayer[:, :, nir_band] - thelayer[:, :, rededge_band]) / (\n", - " thelayer[:, :, nir_band] + thelayer[:, :, rededge_band]\n", - ")\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(np.percentile(masked_ndre, 0.5))\n", - "ndre_hist_max = np.max(np.percentile(masked_ndre, 99.5))\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(\n", - " gamma_corr_rgb,\n", - " masked_ndre,\n", - " figsize=(14, 7),\n", - " title=\"NDRE filtered to only plants over RGB base layer\",\n", - " vmin=min_display_ndre,\n", - " vmax=max_display_ndre,\n", - ")\n", - "fig.savefig(thecapture.uuid + \"_ndre_over_rgb.png\")" - ], - "execution_count": null, - "outputs": [], - "id": "b3fe7a83" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Thermal Imagery\n", - "If our image is from an Altum or Altum-PT 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" - ], - "id": "898cc28a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "if len(thecapture.lw_indices()) > 0:\n", - " lwir_band = thecapture.band_names_lower().index(\"lwir\")\n", - "\n", - " # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n", - " if panchroCam:\n", - " masked_thermal = sharpened_stack[:, :, lwir_band]\n", - " else:\n", - " masked_thermal = im_aligned[:, :, lwir_band]\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", - " # Compute a histogram\n", - " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", - " axis.hist(\n", - " masked_thermal.ravel(),\n", - " bins=512,\n", - " range=(\n", - " np.min(np.percentile(masked_thermal, 1)),\n", - " np.max(np.percentile(masked_thermal, 99)),\n", - " ),\n", - " )\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(\n", - " gamma_corr_rgb,\n", - " masked_thermal,\n", - " figsize=(14, 7),\n", - " title=\"Temperature over True Color\",\n", - " vmin=min_display_therm,\n", - " 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=0.4,\n", - " contour_fmt=\"%.0fC\",\n", - " )\n", - " fig.savefig(thecapture.uuid + \"_thermal_over_rgb.png\")" - ], - "execution_count": null, - "outputs": [], - "id": "176e2c34" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Red vs NIR Reflectance\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." - ], - "id": "bc9b29ec" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "x_band = red_band\n", - "y_band = nir_band\n", - "x_max = np.max(np.percentile(im_aligned[:, :, x_band], 99.99))\n", - "y_max = np.max(np.percentile(im_aligned[:, :, y_band], 99.99))\n", - "\n", - "fig = plt.figure(figsize=(12, 12))\n", - "plt.hexbin(\n", - " im_aligned[:, :, x_band],\n", - " im_aligned[:, :, y_band],\n", - " gridsize=640,\n", - " bins=\"log\",\n", - " extent=(0, x_max, 0, y_max),\n", - ")\n", - "ax = fig.gca()\n", - "ax.set_xlim([0, x_max])\n", - "ax.set_ylim([0, y_max])\n", - "plt.xlabel(\"{} Reflectance\".format(thecapture.band_names()[x_band]))\n", - "plt.ylabel(\"{} Reflectance\".format(thecapture.band_names()[y_band]))\n", - "plt.show()\n", - "\n", - "if panchroCam:\n", - " x_band = red_band\n", - " y_band = nir_band\n", - " x_max = np.max(np.percentile(sharpened_stack[:, :, x_band], 99.99))\n", - " y_max = np.max(np.percentile(sharpened_stack[:, :, y_band], 99.99))\n", - "\n", - " fig = plt.figure(figsize=(12, 12))\n", - " plt.hexbin(\n", - " sharpened_stack[:, :, x_band],\n", - " sharpened_stack[:, :, y_band],\n", - " gridsize=640,\n", - " bins=\"log\",\n", - " extent=(0, x_max, 0, y_max),\n", - " )\n", - " ax = fig.gca()\n", - " ax.set_xlim([0, x_max])\n", - " ax.set_ylim([0, y_max])\n", - " plt.xlabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[x_band]))\n", - " plt.ylabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[y_band]))\n", - " plt.show()" - ], - "execution_count": null, - "outputs": [], - "id": "c0250872" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Print warp_matrices for usage elsewhere, such as Batch Processing\n", - "Lastly, 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." - ], - "id": "4ec2d106" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "if panchroCam:\n", - " print(warp_matrices_SIFT)\n", - "else:\n", - " print(warp_matrices)" - ], - "execution_count": null, - "outputs": [], - "id": "c7cf8f48" - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.12" - } + "cells": [ + { + "cell_type": "markdown", + "id": "573422ea", + "metadata": {}, + "source": [ + "# Active Image Alignment\n", + "For 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 RedEdge (3,M,MX) versions above 3.2.0 which include the \"RigRelatives\" tags, and all RedEdge-P/Altum/Altum-PT imagery. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.\n", + "\n", + "# Opening Images\n", + "As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize all the bands of a MicaSense 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 the `autoreload` extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time." + ] }, - "nbformat": 4, - "nbformat_minor": 5 -} \ No newline at end of file + { + "cell_type": "code", + "execution_count": null, + "id": "b09cd3b1", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7426068b", + "metadata": {}, + "outputs": [], + "source": [ + "import micasense.capture as capture\n", + "\n", + "%matplotlib inline\n", + "from pathlib import Path\n", + "import matplotlib.pyplot as plt\n", + "\n", + "plt.rcParams[\"figure.facecolor\"] = \"w\"\n", + "\n", + "panelNames = None\n", + "\n", + "# set your image paths here. See more here: https://docs.python.org/3/library/pathlib.html\n", + "# if using Windows, you need to an an \"r\" to the path like this: Path(r\"C:\\Files\")\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-MX-DUAL\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0431_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-MX\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0020_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0001_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/ALTUM\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0021_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "# imagePath = Path(\"./data/REDEDGE-P\")\n", + "\n", + "# # these will return lists of image paths as strings\n", + "# imageNames = list(imagePath.glob('IMG_0011_*.tif'))\n", + "# imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "# panelNames = list(imagePath.glob('IMG_0000_*.tif'))\n", + "# panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "imagePath = Path(\"./data/ALTUM-PT\")\n", + "\n", + "# these will return lists of image paths as strings\n", + "imageNames = list(imagePath.glob(\"IMG_0010_*.tif\"))\n", + "imageNames = [x.as_posix() for x in imageNames]\n", + "\n", + "panelNames = list(imagePath.glob(\"IMG_0000_*.tif\"))\n", + "panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "\n", + "if panelNames is not None:\n", + " panelCap = capture.Capture.from_filelist(panelNames)\n", + "else:\n", + " panelCap = None\n", + "\n", + "thecapture = capture.Capture.from_filelist(imageNames)\n", + "\n", + "# get camera model for future use\n", + "cam_model = thecapture.camera_model\n", + "# if this is a multicamera system like the RedEdge-MX Dual,\n", + "# we can combine the two serial numbers to help identify\n", + "# this camera system later.\n", + "if len(thecapture.camera_serials) > 1:\n", + " cam_serial = \"_\".join(thecapture.camera_serials)\n", + " print(cam_serial)\n", + "else:\n", + " cam_serial = thecapture.camera_serial\n", + "\n", + "print(\"Camera model:\", cam_model)\n", + "print(\"Bit depth:\", thecapture.bits_per_pixel)\n", + "print(\"Camera serial number:\", cam_serial)\n", + "print(\"Capture ID:\", thecapture.uuid)\n", + "\n", + "# determine if this sensor has a panchromatic band\n", + "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", + " panchroCam = True\n", + "else:\n", + " panchroCam = False\n", + " panSharpen = False\n", + "\n", + "if panelCap is not None:\n", + " if panelCap.panel_albedo() is not None:\n", + " panel_reflectance_by_band = panelCap.panel_albedo()\n", + " else:\n", + " panel_reflectance_by_band = [0.49] * len(\n", + " thecapture.eo_band_names()\n", + " ) # RedEdge band_index order\n", + " panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)\n", + " irradiance_list = panelCap.panel_irradiance(panel_reflectance_by_band) + [\n", + " 0\n", + " ] # add to account for uncalibrated LWIR band, if applicable\n", + " img_type = \"reflectance\"\n", + " thecapture.plot_undistorted_reflectance(panel_irradiance)\n", + "else:\n", + " if thecapture.dls_present():\n", + " img_type = \"reflectance\"\n", + " irradiance_list = thecapture.dls_irradiance() + [0]\n", + " thecapture.plot_undistorted_reflectance(thecapture.dls_irradiance())\n", + " else:\n", + " img_type = \"radiance\"\n", + " thecapture.plot_undistorted_radiance()\n", + " irradiance_list = None" + ] + }, + { + "cell_type": "markdown", + "id": "8de5a774", + "metadata": {}, + "source": [ + "# Check for existing warp matrices \n", + "If we have already successfully aligned captures from this specific camera, we can typically save some time and use the alignment warp matrices for other captures from the same camera" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35de1be0", + "metadata": {}, + "outputs": [], + "source": [ + "from micasense.warp_io import load_warp_matrices\n", + "\n", + "if panchroCam:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", + "else:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", + "\n", + "if Path(\"./\" + warp_matrices_filename).is_file():\n", + " print(\"Found existing warp matrices for camera\", cam_serial)\n", + " loaded = load_warp_matrices(warp_matrices_filename, as_projective=panchroCam)\n", + " print(\"Warp matrices successfully loaded.\")\n", + "\n", + " if panchroCam:\n", + " warp_matrices_SIFT = loaded\n", + " else:\n", + " warp_matrices = loaded\n", + "else:\n", + " print(\"No existing warp matrices found. Create them later in the notebook.\")\n", + " warp_matrices_SIFT = False\n", + " warp_matrices = False" + ] + }, + { + "cell_type": "markdown", + "id": "185627ff", + "metadata": {}, + "source": [ + "# Unwarp and Align (OpenCV method for RedEdge3/M/MX/Dual and original Altum)\n", + "Alignment is a three-step process:\n", + "\n", + "Images are unwarped using the built-in lens calibration\n", + "A transformation is found to align each band to a common band\n", + "The aligned images are combined and cropped, removing pixels which don't overlap in all bands.\n", + "We provide utilities to find the alignment transformations within a single capture. Our experience shows that once a good alignment 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 alignment 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, + "id": "e154cdd1", + "metadata": {}, + "outputs": [], + "source": [ + "import cv2\n", + "import time\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", + "# We use a different alignment method for RedEdge-P and Altum-PT,\n", + "# so if the imagery is from this kind of camera, skip this step\n", + "if not panchroCam:\n", + " st = time.time()\n", + " # set to True if you'd like to ignore existing warp matrices and create new ones\n", + " regenerate = True\n", + " pyramid_levels = (\n", + " 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment\n", + " )\n", + " max_alignment_iterations = 10\n", + "\n", + " # match_index:\n", + " # for non-panchromatic cameras we want to use band 1, which is green\n", + " # the green band has zero rig relative offsets\n", + " # NOTE: These band numbers are zero-indexed\n", + " # special parameters for RedEdge-MX dual camera system\n", + " if len(thecapture.eo_band_names()) == 10:\n", + " print(\"is 10 band\")\n", + " pyramid_levels = 3\n", + " match_index = 4\n", + " max_alignment_iterations = 20\n", + " else:\n", + " match_index = 1\n", + "\n", + " warp_mode = (\n", + " cv2.MOTION_HOMOGRAPHY\n", + " ) # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY\n", + "\n", + " print(\n", + " \"Aligning images. Depending on settings this can take from a few seconds to many minutes\"\n", + " )\n", + " # Can potentially increase max_iterations for better results, but longer runtimes\n", + " if warp_matrices and not regenerate:\n", + " print(\"Using existing warp matrices...\")\n", + " try:\n", + " irradiance = panel_irradiance + [0]\n", + " except NameError:\n", + " irradiance = None\n", + " else:\n", + " warp_matrices, alignment_pairs = imageutils.align_capture(\n", + " thecapture,\n", + " ref_index=match_index,\n", + " max_iterations=max_alignment_iterations,\n", + " warp_mode=warp_mode,\n", + " pyramid_levels=pyramid_levels,\n", + " )\n", + "\n", + " print(\"Finished Aligning\")\n", + " et = time.time()\n", + " elapsed_time = et - st\n", + " print(\"Alignment time:\", int(elapsed_time), \"seconds\")" + ] + }, + { + "cell_type": "markdown", + "id": "c06bfdf3", + "metadata": {}, + "source": [ + "# Crop Aligned Images and create aligned capture stack \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 take 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 to 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, + "id": "5fbfb346", + "metadata": {}, + "outputs": [], + "source": [ + "if not panchroCam:\n", + " cropped_dimensions, edges = imageutils.find_crop_bounds(\n", + " thecapture, warp_matrices, warp_mode=warp_mode, reference_band=match_index\n", + " )\n", + " print(\"Cropped dimensions:\", cropped_dimensions)\n", + " im_aligned = thecapture.create_aligned_capture(\n", + " warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "c362623a", + "metadata": {}, + "source": [ + "# Alignment and pan-sharpening for RedEdge-P and Altum-PT\n", + "For older cameras we use OpenCV for capture alignment. For RedEdge-P and Altum-PT, we use SIFT (scale-invariant feature transform). We will use SIFT to create the warp matrices, then use these warp matrices to align the capture. See more here: \n", + "https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_sift.html\n", + "https://en.wikipedia.org/wiki/Scale-invariant_feature_transform\n", + "\n", + "For sensors with a panchromatic band (Altum-PT or RedEdge-P), we may wish to create a pan-sharpened stack. This example uses a linear interpolation method for pan-sharpening.\n", + "\n", + "The `radiometric_pan_sharpen` function takes in the SIFT warp matrices and outputs an upsampled image stack (all bands are changed to the resolution of the panchromatic sensor) and a pan-sharpened stack. \n", + "\n", + "Note: it is important to choose an initial alignment image that has a lot of straight, manmade features, such as roads or buildings. Trying to align a capture that is strictly vegetation will take a long time and may not be very good. This process can take a while, depending on how feature-rich the capture is. We have seen some Altum-PT captures align after 2 minutes, and some can take upwards of 30 minutes. Be patient! Once an alignment has been found, it can be used on other captures from the same camera, and the alignment process will be much faster. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12c2a278", + "metadata": {}, + "outputs": [], + "source": [ + "# from micasense.imageutils import brovey_pan_sharpen,radiometric_pan_sharpen\n", + "from skimage.transform import ProjectiveTransform\n", + "import time\n", + "\n", + "if panchroCam:\n", + " # set to True if you'd like to ignore existing warp matrices and create new ones\n", + " regenerate = True\n", + " st = time.time()\n", + " if not warp_matrices_SIFT or regenerate:\n", + " print(\"Generating new warp matrices...\")\n", + " warp_matrices_SIFT = thecapture.SIFT_align_capture(min_matches=10)\n", + "\n", + " sharpened_stack, upsampled = thecapture.radiometric_pan_sharpened_aligned_capture(\n", + " warp_matrices=warp_matrices_SIFT,\n", + " irradiance_list=irradiance_list,\n", + " img_type=img_type,\n", + " )\n", + "\n", + " # we can also use the Rig Relatives from the image metadata to do a quick, rudimentary alignment\n", + " # warp_matrices0=thecapture.get_warp_matrices(ref_index=5)\n", + " # sharpened_stack,upsampled = radiometric_pan_sharpen(thecapture,warp_matrices=warp_matrices0)\n", + "\n", + " print(\"Pansharpened shape:\", sharpened_stack.shape)\n", + " print(\"Upsampled shape:\", upsampled.shape)\n", + " # re-assign to im_aligned to match rest of code\n", + " im_aligned = upsampled\n", + " et = time.time()\n", + " elapsed_time = et - st\n", + " print(\"Alignment and pan-sharpening time:\", int(elapsed_time), \"seconds\")" + ] + }, + { + "cell_type": "markdown", + "id": "3e9d51b1", + "metadata": {}, + "source": [ + "# Save warp matrices\n", + "Once an alignment for your camera has been found, it can be saved to a file for later use with this notebook. It can also be used on the Batch Alignment notebook for aligning all of the captures from an entire flight." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84855db4", + "metadata": {}, + "outputs": [], + "source": [ + "from micasense.warp_io import save_warp_matrices\n", + "\n", + "if panchroCam:\n", + " working_wm = warp_matrices_SIFT\n", + "else:\n", + " working_wm = warp_matrices\n", + "if not Path(\"./\" + warp_matrices_filename).is_file() or regenerate:\n", + " save_warp_matrices(warp_matrices_filename, working_wm)\n", + " print(\"Saved to\", Path(\"./\" + warp_matrices_filename).resolve())\n", + "else:\n", + " print(\"Matrices already exist at\", Path(\"./\" + warp_matrices_filename).resolve())" + ] + }, + { + "cell_type": "markdown", + "id": "47d3fb77", + "metadata": {}, + "source": [ + "# Multispectral band histogram \n", + "We can compare the radiance between multispectral bands, and between upsampled/pansharpened in the case of RedEdge-P and Altum-PT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e036b45", + "metadata": {}, + "outputs": [], + "source": [ + "theColors = {\n", + " \"Blue\": \"blue\",\n", + " \"Green\": \"green\",\n", + " \"Red\": \"red\",\n", + " \"Red edge\": \"maroon\",\n", + " \"NIR\": \"purple\",\n", + " \"Panchro\": \"yellow\",\n", + " \"PanchroB\": \"orange\",\n", + " \"Red edge-740\": \"salmon\",\n", + " \"Red Edge\": \"maroon\",\n", + " \"Blue-444\": \"aqua\",\n", + " \"Green-531\": \"lime\",\n", + " \"Red-650\": \"lightcoral\",\n", + " \"Red edge-705\": \"brown\",\n", + "}\n", + "\n", + "eo_count = len(thecapture.eo_indices())\n", + "multispec_min = np.min(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 0.01))\n", + "multispec_max = np.max(np.percentile(im_aligned[:, :, 1:eo_count].flatten(), 99.99))\n", + "\n", + "theRange = (multispec_min, multispec_max)\n", + "\n", + "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + "for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", + " axis.hist(\n", + " im_aligned[:, :, x].ravel(),\n", + " bins=512,\n", + " range=theRange,\n", + " histtype=\"step\",\n", + " label=y,\n", + " color=theColors[y],\n", + " linewidth=1.5,\n", + " )\n", + "plt.title(\"Multispectral histogram (radiance)\")\n", + "axis.legend()\n", + "plt.show()\n", + "\n", + "if panchroCam:\n", + " eo_count = len(thecapture.eo_indices())\n", + " multispec_min = np.min(\n", + " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 0.01)\n", + " )\n", + " multispec_max = np.max(\n", + " np.percentile(sharpened_stack[:, :, 1:eo_count].flatten(), 99.99)\n", + " )\n", + "\n", + " theRange = (multispec_min, multispec_max)\n", + "\n", + " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + " for x, y in zip(thecapture.eo_indices(), thecapture.eo_band_names()):\n", + " axis.hist(\n", + " sharpened_stack[:, :, x].ravel(),\n", + " bins=512,\n", + " range=theRange,\n", + " histtype=\"step\",\n", + " label=y,\n", + " color=theColors[y],\n", + " linewidth=1.5,\n", + " )\n", + " plt.title(\"Pan-sharpened multispectral histogram (radiance)\")\n", + " axis.legend()\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "3ccb819b", + "metadata": {}, + "source": [ + "# Visualize Aligned Images\n", + "Once the transformation has been found, it can be verified by compositing 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 external software such as QGIS. Useful composites are a naturally colored RGB as well as color infrared, or CIR." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8626c90b", + "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 = [\n", + " thecapture.band_names_lower().index(\"red\"),\n", + " thecapture.band_names_lower().index(\"green\"),\n", + " thecapture.band_names_lower().index(\"blue\"),\n", + "]\n", + "cir_band_indices = [\n", + " thecapture.band_names_lower().index(\"nir\"),\n", + " thecapture.band_names_lower().index(\"red\"),\n", + " thecapture.band_names_lower().index(\"green\"),\n", + "]\n", + "\n", + "# Create normalized stacks for viewing\n", + "im_display = np.zeros(\n", + " (im_aligned.shape[0], im_aligned.shape[1], im_aligned.shape[2]), dtype=float\n", + ")\n", + "im_min = np.percentile(\n", + " im_aligned[:, :, rgb_band_indices].flatten(), 0.5\n", + ") # modify these percentiles to adjust contrast\n", + "im_max = np.percentile(\n", + " im_aligned[:, :, rgb_band_indices].flatten(), 99.5\n", + ") # for many images, 0.5 and 99.5 are good values\n", + "\n", + "if panchroCam:\n", + " im_display_sharp = np.zeros(\n", + " (sharpened_stack.shape[0], sharpened_stack.shape[1], sharpened_stack.shape[2]),\n", + " dtype=float,\n", + " )\n", + " im_min_sharp = np.percentile(\n", + " sharpened_stack[:, :, rgb_band_indices].flatten(), 0.5\n", + " ) # modify these percentiles to adjust contrast\n", + " im_max_sharp = np.percentile(\n", + " sharpened_stack[:, :, rgb_band_indices].flatten(), 99.5\n", + " ) # for many images, 0.5 and 99.5 are good values\n", + "\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", + " if panchroCam:\n", + " im_display_sharp[:, :, i] = imageutils.normalize(\n", + " sharpened_stack[:, :, i], im_min_sharp, im_max_sharp\n", + " )\n", + "\n", + "rgb = im_display[:, :, rgb_band_indices]\n", + "\n", + "if panchroCam:\n", + " rgb_sharp = im_display_sharp[:, :, rgb_band_indices]\n", + "\n", + "nir_band = thecapture.band_names_lower().index(\"nir\")\n", + "red_band = thecapture.band_names_lower().index(\"red\")\n", + "\n", + "ndvi = (im_aligned[:, :, nir_band] - im_aligned[:, :, red_band]) / (\n", + " im_aligned[:, :, nir_band] + im_aligned[:, :, red_band]\n", + ")\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", + "if panchroCam:\n", + " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)\n", + "else:\n", + " fig, ax1 = plt.subplots(1, 1, figsize=figsize)\n", + "ax1.set_title(\"Red-Green-Blue Composite\")\n", + "ax1.imshow(rgb)\n", + "if panchroCam:\n", + " ax2.set_title(\"Red-Green-Blue Composite (pan-sharpened)\")\n", + " ax2.imshow(rgb_sharp)\n", + "\n", + "fig, (ax3, ax4) = plt.subplots(1, 2, figsize=figsize)\n", + "ax3.set_title(\"NDVI\")\n", + "ax3.imshow(ndvi)\n", + "ax4.set_title(\"Color Infrared (CIR) Composite\")\n", + "ax4.imshow(cir)\n", + "\n", + "# set custom lims if you want to zoom in to image to see more detail\n", + "# this is useful for comparing upsampled and pan-sharpened stacks\n", + "# custom_xlim=(1500,2000)\n", + "# custom_ylim=(2000,1500)\n", + "# plt.setp([ax1,ax2,ax3,ax4], xlim=custom_xlim, ylim=custom_ylim)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "17f71519", + "metadata": {}, + "source": [ + "# Image Enhancement\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, + "id": "1185fd16", + "metadata": {}, + "outputs": [], + "source": [ + "if panchroCam:\n", + " rgb = rgb_sharp\n", + "# 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", + "id": "6d1f6189", + "metadata": {}, + "source": [ + "# Stack Export\n", + "We can easily export the stacks into an image 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. The stacks include geographic information. \n", + "\n", + "If you prefer, you may set `sort_by_wavelength` to `True` in the `save_capture_as_stack` function.\n", + "\n", + "Unless otherwise specified, this will save in your working folder, that is your `imageprocessing` directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2d716be", + "metadata": {}, + "outputs": [], + "source": [ + "# set output name to unique capture ID, e.g. FWoNSvgDNBX63Xv378qs\n", + "outputName = thecapture.uuid\n", + "\n", + "st = time.time()\n", + "if panchroCam:\n", + " # in this example, we can export both a pan-sharpened stack and an upsampled stack\n", + " # so you can compare them in GIS. In practice, you would typically only output the pansharpened stack\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-pansharpened.tif\", sort_by_wavelength=True, pansharpen=True\n", + " )\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-upsampled.tif\", sort_by_wavelength=True, pansharpen=False\n", + " )\n", + "else:\n", + " thecapture.save_capture_as_stack(\n", + " outputName + \"-noPanels.tif\", sort_by_wavelength=True\n", + " )\n", + "\n", + "et = time.time()\n", + "elapsed_time = et - st\n", + "print(\"Time to save stacks:\", int(elapsed_time), \"seconds.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f4d9cd64", + "metadata": {}, + "source": [ + "# NDVI Computation\n", + "For raw index computation on single images, the `numpy` package provides a simple way to do math and simple visualization on images. Below, we compute and visualize an image histogram, and then use that to pick a color map range for visualizing the NDVI of an image.\n", + "\n", + "## Plant Classification\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, + "id": "4a572a55", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "nir_band = thecapture.band_names_lower().index(\"nir\")\n", + "red_band = thecapture.band_names_lower().index(\"red\")\n", + "\n", + "thelayer = im_aligned\n", + "if panchroCam:\n", + " thelayer = sharpened_stack\n", + "np.seterr(\n", + " divide=\"ignore\", invalid=\"ignore\"\n", + ") # 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 = (thelayer[:, :, nir_band] - thelayer[:, :, red_band]) / (\n", + " thelayer[:, :, nir_band] + thelayer[:, :, red_band]\n", + ")\n", + "print(\"Image type:\", img_type)\n", + "\n", + "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n", + "# this does not seem to work on panchro stacks\n", + "if img_type == \"reflectance\":\n", + " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < 0.20, ndvi)\n", + "elif img_type == \"radiance\":\n", + " lower_pct_radiance = np.percentile(thelayer[:, :, nir_band], 10.0)\n", + " ndvi = np.ma.masked_where(thelayer[:, :, nir_band] < lower_pct_radiance, ndvi)\n", + "# Compute and display a histogram\n", + "# ndvi_hist_min = np.min(ndvi)\n", + "# ndvi_hist_max = np.max(ndvi)\n", + "ndvi_hist_min = np.min(np.percentile(ndvi, 0.5))\n", + "ndvi_hist_max = np.max(np.percentile(ndvi, 99.5))\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(\n", + " ndvi.flatten(), 99.5\n", + ") # 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(\n", + " gamma_corr_rgb,\n", + " masked_ndvi,\n", + " figsize=(14, 7),\n", + " title=\"NDVI filtered to only plants over RGB base layer\",\n", + " vmin=min_display_ndvi,\n", + " vmax=max_display_ndvi,\n", + ")\n", + "fig.savefig(thecapture.uuid + \"_ndvi_over_rgb.png\")" + ] + }, + { + "cell_type": "markdown", + "id": "d8b4cf19", + "metadata": {}, + "source": [ + "# NDRE Computation\n", + "In the same manner, we can compute, filter, and display another index useful for MicaSense cameras, 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, + "id": "b3fe7a83", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n", + "rededge_band = thecapture.band_names_lower().index(\"red edge\")\n", + "ndre = (thelayer[:, :, nir_band] - thelayer[:, :, rededge_band]) / (\n", + " thelayer[:, :, nir_band] + thelayer[:, :, rededge_band]\n", + ")\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(np.percentile(masked_ndre, 0.5))\n", + "ndre_hist_max = np.max(np.percentile(masked_ndre, 99.5))\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(\n", + " gamma_corr_rgb,\n", + " masked_ndre,\n", + " figsize=(14, 7),\n", + " title=\"NDRE filtered to only plants over RGB base layer\",\n", + " vmin=min_display_ndre,\n", + " vmax=max_display_ndre,\n", + ")\n", + "fig.savefig(thecapture.uuid + \"_ndre_over_rgb.png\")" + ] + }, + { + "cell_type": "markdown", + "id": "898cc28a", + "metadata": {}, + "source": [ + "# Thermal Imagery\n", + "If our image is from an Altum or Altum-PT 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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "176e2c34", + "metadata": {}, + "outputs": [], + "source": [ + "if len(thecapture.lw_indices()) > 0:\n", + " lwir_band = thecapture.band_names_lower().index(\"lwir\")\n", + "\n", + " # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n", + " if panchroCam:\n", + " masked_thermal = sharpened_stack[:, :, lwir_band]\n", + " else:\n", + " masked_thermal = im_aligned[:, :, lwir_band]\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", + " # Compute a histogram\n", + " fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + " axis.hist(\n", + " masked_thermal.ravel(),\n", + " bins=512,\n", + " range=(\n", + " np.min(np.percentile(masked_thermal, 1)),\n", + " np.max(np.percentile(masked_thermal, 99)),\n", + " ),\n", + " )\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(\n", + " gamma_corr_rgb,\n", + " masked_thermal,\n", + " figsize=(14, 7),\n", + " title=\"Temperature over True Color\",\n", + " vmin=min_display_therm,\n", + " 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=0.4,\n", + " contour_fmt=\"%.0fC\",\n", + " )\n", + " fig.savefig(thecapture.uuid + \"_thermal_over_rgb.png\")" + ] + }, + { + "cell_type": "markdown", + "id": "bc9b29ec", + "metadata": {}, + "source": [ + "# Red vs NIR Reflectance\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, + "id": "c0250872", + "metadata": {}, + "outputs": [], + "source": [ + "x_band = red_band\n", + "y_band = nir_band\n", + "x_max = np.max(np.percentile(im_aligned[:, :, x_band], 99.99))\n", + "y_max = np.max(np.percentile(im_aligned[:, :, y_band], 99.99))\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "plt.hexbin(\n", + " im_aligned[:, :, x_band],\n", + " im_aligned[:, :, y_band],\n", + " gridsize=640,\n", + " bins=\"log\",\n", + " extent=(0, x_max, 0, y_max),\n", + ")\n", + "ax = fig.gca()\n", + "ax.set_xlim([0, x_max])\n", + "ax.set_ylim([0, y_max])\n", + "plt.xlabel(\"{} Reflectance\".format(thecapture.band_names()[x_band]))\n", + "plt.ylabel(\"{} Reflectance\".format(thecapture.band_names()[y_band]))\n", + "plt.show()\n", + "\n", + "if panchroCam:\n", + " x_band = red_band\n", + " y_band = nir_band\n", + " x_max = np.max(np.percentile(sharpened_stack[:, :, x_band], 99.99))\n", + " y_max = np.max(np.percentile(sharpened_stack[:, :, y_band], 99.99))\n", + "\n", + " fig = plt.figure(figsize=(12, 12))\n", + " plt.hexbin(\n", + " sharpened_stack[:, :, x_band],\n", + " sharpened_stack[:, :, y_band],\n", + " gridsize=640,\n", + " bins=\"log\",\n", + " extent=(0, x_max, 0, y_max),\n", + " )\n", + " ax = fig.gca()\n", + " ax.set_xlim([0, x_max])\n", + " ax.set_ylim([0, y_max])\n", + " plt.xlabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[x_band]))\n", + " plt.ylabel(\"{} Reflectance (pan-sharpened)\".format(thecapture.band_names()[y_band]))\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "4ec2d106", + "metadata": {}, + "source": [ + "# Print warp_matrices for usage elsewhere, such as Batch Processing\n", + "Lastly, 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, + "id": "c7cf8f48", + "metadata": {}, + "outputs": [], + "source": [ + "if panchroCam:\n", + " print(warp_matrices_SIFT)\n", + "else:\n", + " print(warp_matrices)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Batch Processing v2.ipynb b/Batch Processing v2.ipynb index 969934f..db6f951 100644 --- a/Batch Processing v2.ipynb +++ b/Batch Processing v2.ipynb @@ -1,379 +1,378 @@ { - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Batch Processing Example\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. Part of this process (via `imageutils.write_exif_to_stack`) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.\n", - "\n", - "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support" - ], - "id": "20e37d9b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ], - "execution_count": null, - "outputs": [], - "id": "c64ead96" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Load Images into ImageSet\n" - ], - "id": "80646655" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "from ipywidgets import FloatProgress, Layout\n", - "from IPython.display import display\n", - "import micasense.capture as capture\n", - "import os\n", - "from pathlib import Path\n", - "\n", - "# set to True if you have an Altum-PT\n", - "# or RedEdge-P and wish to output pan-sharpened stacks\n", - "panSharpen = True\n", - "\n", - "# If creating a lot of stacks, it is more efficient to save the metadata\n", - "# and then write all of the exif to the images after the stacks are created\n", - "write_exif_to_individual_stacks = False\n", - "\n", - "panelNames = None\n", - "useDLS = True\n", - "\n", - "# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html\n", - "imagePath = Path(\"./data/REDEDGE-MX\")\n", - "\n", - "# these will return lists of image paths as strings. Comment out of you aren't using panels.\n", - "panelNames = list(imagePath.glob(\"IMG_0001_*.tif\"))\n", - "panelNames = [x.as_posix() for x in panelNames]\n", - "\n", - "if panelNames:\n", - " panelCap = capture.Capture.from_filelist(panelNames)\n", - "\n", - "# destinations on your computer to put the stacks\n", - "# and RGB thumbnails\n", - "outputPath = imagePath / \"..\" / \"stacks\"\n", - "thumbnailPath = outputPath / \"thumbnails\"\n", - "\n", - "cam_model = panelCap.camera_model\n", - "cam_serial = panelCap.camera_serial\n", - "\n", - "# determine if this sensor has a panchromatic band\n", - "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", - " panchroCam = True\n", - "else:\n", - " panchroCam = False\n", - " panSharpen = False\n", - "\n", - "# if this is a multicamera system like the RedEdge-MX Dual,\n", - "# we can combine the two serial numbers to help identify\n", - "# this camera system later.\n", - "if len(panelCap.camera_serials) > 1:\n", - " cam_serial = \"_\".join(panelCap.camera_serials)\n", - " print(\"Serial number:\", cam_serial)\n", - "else:\n", - " cam_serial = panelCap.camera_serial\n", - " print(\"Serial number:\", cam_serial)\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; but 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(\n", - " v is None for v in panelCap.panel_albedo()\n", - " ):\n", - " panel_reflectance_by_band = panelCap.panel_albedo()\n", - " else:\n", - " panel_reflectance_by_band = [0.49] * len(\n", - " panelCap.eo_band_names()\n", - " ) # 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\"" - ], - "execution_count": null, - "outputs": [], - "id": "2081b06b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import micasense.imageset as imageset\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", - "\n", - "\n", - "def update_f(val):\n", - " if (\n", - " val - f.value\n", - " ) > 0.005 or val == 1: # reduces cpu usage from updating the progressbar by 10x\n", - " f.value = val\n", - "\n", - "# time is not recognized and would remove the imageset as unused\n", - "imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", - "update_f(1.0)" - ], - "execution_count": null, - "outputs": [], - "id": "60736009" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Capture map\n", - "We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath." - ], - "id": "6c0a234d" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import numpy as np\n", - "from mapboxgl.viz import CircleViz\n", - "from mapboxgl.utils import df_to_geojson\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 access token here (https://account.mapbox.com/access-tokens/)\n", - "token = \"YOUR_MAPBOX_ACCESS_TOKEN\"\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", - "\n", - "geojson_data = df_to_geojson(df, columns[3:], lat=\"latitude\", lon=\"longitude\")\n", - "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", - "color_stops = create_color_stops(breaks, colors=\"YlOrRd\")\n", - "\n", - "viz = CircleViz(\n", - " geojson_data,\n", - " access_token=token,\n", - " color_property=color_property,\n", - " color_stops=color_stops,\n", - " center=[df[\"longitude\"].median(), df[\"latitude\"].median()],\n", - " zoom=16,\n", - " height=\"600px\",\n", - " style=\"mapbox://styles/mapbox/satellite-streets-v9\",\n", - ")\n", - "viz.show()" - ], - "execution_count": null, - "outputs": [], - "id": "9e9c437c" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Define which warp method to use\n", - "For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) 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. To use this simple alignment, simply set `warp_matrices=None` \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." - ], - "id": "e540a655" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "from micasense.warp_io import load_warp_matrices\n", - "\n", - "if panchroCam:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", - "else:\n", - " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", - "\n", - "if Path(\"./\" + warp_matrices_filename).is_file():\n", - " print(\"Found existing warp matrices for camera\", cam_serial)\n", - " loaded = load_warp_matrices(\n", - " warp_matrices_filename, as_projective=panchroCam\n", - " )\n", - "\n", - " if panchroCam:\n", - " warp_matrices_SIFT = loaded\n", - " else:\n", - " warp_matrices = loaded\n", - " print(\"Loaded warp matrices from\", Path(\"./\" + warp_matrices_filename).resolve())\n", - "else:\n", - " print(\"No warp matrices found at expected location:\", warp_matrices_filename)" - ], - "execution_count": null, - "outputs": [], - "id": "635ad8d8" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Align images and save each capture to a layered TIFF file" - ], - "id": "3aae7f09" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import datetime\n", - "import micasense.imageutils as imageutils\n", - "\n", - "exif_list = []\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", - "\n", - "\n", - "def update_f2(val):\n", - " f2.value = val\n", - "\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 = str(i).zfill(4) + \"_\" + capture.uuid + \".tif\"\n", - " thumbnailFilename = str(i).zfill(4) + \"_\" + 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", - " if panchroCam:\n", - " capture.radiometric_pan_sharpened_aligned_capture(\n", - " warp_matrices=warp_matrices_SIFT,\n", - " irradiance_list=capture.dls_irradiace(),\n", - " img_type=img_type,\n", - " write_exif=write_exif_to_individual_stacks,\n", - " )\n", - " else:\n", - " capture.create_aligned_capture(\n", - " irradiance_list=irradiance, warp_matrices=warp_matrices\n", - " )\n", - " exif_list.append(\n", - " imageutils.prepare_exif_for_stacks(capture, fullOutputPath)\n", - " )\n", - " capture.save_capture_as_stack(\n", - " fullOutputPath,\n", - " pansharpen=panSharpen,\n", - " sort_by_wavelength=True,\n", - " write_exif=write_exif_to_individual_stacks,\n", - " )\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(\n", - " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", - " float(len(imgset.captures)) / float((end - start).total_seconds())\n", - " )\n", - ")" - ], - "execution_count": null, - "outputs": [], - "id": "cdc195da" - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Write EXIF data to stacks\n", - "As mentioned above, it is more time intensive to write the exif data to each image as it is created. Here, we write the exif data after all of the TIFF files have been created. This should take a few seconds per stack." - ], - "id": "f4b88c80" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "if not write_exif_to_individual_stacks:\n", - " start = datetime.datetime.now()\n", - " for exif in exif_list:\n", - " imageutils.write_exif_to_stack(existing_exif_list=exif)\n", - " end = datetime.datetime.now()\n", - " print(\"Saving time: {}\".format(end - start))\n", - " print(\n", - " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", - " float(len(exif_list)) / float((end - start).total_seconds())\n", - " )\n", - " )" - ], - "execution_count": null, - "outputs": [], - "id": "fac88b6a" - } - ], - "metadata": { - "kernelspec": { - "display_name": "imageprocessing", - "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.12.12" - } + "cells": [ + { + "cell_type": "markdown", + "id": "20e37d9b", + "metadata": {}, + "source": [ + "# Batch Processing Example\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. Part of this process (via `imageutils.write_exif_to_stack`) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.\n", + "\n", + "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support" + ] }, - "nbformat": 4, - "nbformat_minor": 5 -} \ No newline at end of file + { + "cell_type": "code", + "execution_count": null, + "id": "c64ead96", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "id": "80646655", + "metadata": {}, + "source": [ + "# Load Images into ImageSet\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2081b06b", + "metadata": {}, + "outputs": [], + "source": [ + "from ipywidgets import FloatProgress, Layout\n", + "from IPython.display import display\n", + "import micasense.capture as capture\n", + "import os\n", + "from pathlib import Path\n", + "\n", + "# set to True if you have an Altum-PT\n", + "# or RedEdge-P and wish to output pan-sharpened stacks\n", + "panSharpen = True\n", + "\n", + "# If creating a lot of stacks, it is more efficient to save the metadata\n", + "# and then write all of the exif to the images after the stacks are created\n", + "write_exif_to_individual_stacks = False\n", + "\n", + "panelNames = None\n", + "useDLS = True\n", + "\n", + "# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html\n", + "imagePath = Path(\"./data/REDEDGE-MX\")\n", + "\n", + "# these will return lists of image paths as strings. Comment out of you aren't using panels.\n", + "panelNames = list(imagePath.glob(\"IMG_0001_*.tif\"))\n", + "panelNames = [x.as_posix() for x in panelNames]\n", + "\n", + "if panelNames:\n", + " panelCap = capture.Capture.from_filelist(panelNames)\n", + "\n", + "# destinations on your computer to put the stacks\n", + "# and RGB thumbnails\n", + "outputPath = imagePath / \"..\" / \"stacks\"\n", + "thumbnailPath = outputPath / \"thumbnails\"\n", + "\n", + "cam_model = panelCap.camera_model\n", + "cam_serial = panelCap.camera_serial\n", + "\n", + "# determine if this sensor has a panchromatic band\n", + "if cam_model == \"RedEdge-P\" or cam_model == \"Altum-PT\":\n", + " panchroCam = True\n", + "else:\n", + " panchroCam = False\n", + " panSharpen = False\n", + "\n", + "# if this is a multicamera system like the RedEdge-MX Dual,\n", + "# we can combine the two serial numbers to help identify\n", + "# this camera system later.\n", + "if len(panelCap.camera_serials) > 1:\n", + " cam_serial = \"_\".join(panelCap.camera_serials)\n", + " print(\"Serial number:\", cam_serial)\n", + "else:\n", + " cam_serial = panelCap.camera_serial\n", + " print(\"Serial number:\", cam_serial)\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; but 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(\n", + " v is None for v in panelCap.panel_albedo()\n", + " ):\n", + " panel_reflectance_by_band = panelCap.panel_albedo()\n", + " else:\n", + " panel_reflectance_by_band = [0.49] * len(\n", + " panelCap.eo_band_names()\n", + " ) # 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, + "id": "60736009", + "metadata": {}, + "outputs": [], + "source": [ + "import micasense.imageset as imageset\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", + "\n", + "\n", + "def update_f(val):\n", + " if (\n", + " val - f.value\n", + " ) > 0.005 or val == 1: # reduces cpu usage from updating the progressbar by 10x\n", + " f.value = val\n", + "\n", + "\n", + "# time is not recognized and would remove the imageset as unused\n", + "imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", + "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n", + "update_f(1.0)" + ] + }, + { + "cell_type": "markdown", + "id": "6c0a234d", + "metadata": {}, + "source": [ + "# Capture map\n", + "We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e9c437c", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from mapboxgl.viz import CircleViz\n", + "from mapboxgl.utils import df_to_geojson\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 access token here (https://account.mapbox.com/access-tokens/)\n", + "token = \"YOUR_MAPBOX_ACCESS_TOKEN\"\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", + "\n", + "geojson_data = df_to_geojson(df, columns[3:], lat=\"latitude\", lon=\"longitude\")\n", + "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n", + "color_stops = create_color_stops(breaks, colors=\"YlOrRd\")\n", + "\n", + "viz = CircleViz(\n", + " geojson_data,\n", + " access_token=token,\n", + " color_property=color_property,\n", + " color_stops=color_stops,\n", + " center=[df[\"longitude\"].median(), df[\"latitude\"].median()],\n", + " zoom=16,\n", + " height=\"600px\",\n", + " style=\"mapbox://styles/mapbox/satellite-streets-v9\",\n", + ")\n", + "viz.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e540a655", + "metadata": {}, + "source": [ + "# Define which warp method to use\n", + "For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) 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. To use this simple alignment, simply set `warp_matrices=None` \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, + "id": "635ad8d8", + "metadata": {}, + "outputs": [], + "source": [ + "from micasense.warp_io import load_warp_matrices\n", + "\n", + "if panchroCam:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n", + "else:\n", + " warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n", + "\n", + "if Path(\"./\" + warp_matrices_filename).is_file():\n", + " print(\"Found existing warp matrices for camera\", cam_serial)\n", + " loaded = load_warp_matrices(warp_matrices_filename, as_projective=panchroCam)\n", + "\n", + " if panchroCam:\n", + " warp_matrices_SIFT = loaded\n", + " else:\n", + " warp_matrices = loaded\n", + " print(\"Loaded warp matrices from\", Path(\"./\" + warp_matrices_filename).resolve())\n", + "else:\n", + " print(\"No warp matrices found at expected location:\", warp_matrices_filename)" + ] + }, + { + "cell_type": "markdown", + "id": "3aae7f09", + "metadata": {}, + "source": [ + "## Align images and save each capture to a layered TIFF file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdc195da", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "import micasense.imageutils as imageutils\n", + "\n", + "exif_list = []\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", + "\n", + "\n", + "def update_f2(val):\n", + " f2.value = val\n", + "\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 = str(i).zfill(4) + \"_\" + capture.uuid + \".tif\"\n", + " thumbnailFilename = str(i).zfill(4) + \"_\" + 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", + " if panchroCam:\n", + " capture.radiometric_pan_sharpened_aligned_capture(\n", + " warp_matrices=warp_matrices_SIFT,\n", + " irradiance_list=capture.dls_irradiace(),\n", + " img_type=img_type,\n", + " write_exif=write_exif_to_individual_stacks,\n", + " )\n", + " else:\n", + " capture.create_aligned_capture(\n", + " irradiance_list=irradiance, warp_matrices=warp_matrices\n", + " )\n", + " exif_list.append(\n", + " imageutils.prepare_exif_for_stacks(capture, fullOutputPath)\n", + " )\n", + " capture.save_capture_as_stack(\n", + " fullOutputPath,\n", + " pansharpen=panSharpen,\n", + " sort_by_wavelength=True,\n", + " write_exif=write_exif_to_individual_stacks,\n", + " )\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(\n", + " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", + " float(len(imgset.captures)) / float((end - start).total_seconds())\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "f4b88c80", + "metadata": {}, + "source": [ + "# Write EXIF data to stacks\n", + "As mentioned above, it is more time intensive to write the exif data to each image as it is created. Here, we write the exif data after all of the TIFF files have been created. This should take a few seconds per stack." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fac88b6a", + "metadata": {}, + "outputs": [], + "source": [ + "if not write_exif_to_individual_stacks:\n", + " start = datetime.datetime.now()\n", + " for exif in exif_list:\n", + " imageutils.write_exif_to_stack(existing_exif_list=exif)\n", + " end = datetime.datetime.now()\n", + " print(\"Saving time: {}\".format(end - start))\n", + " print(\n", + " \"Alignment+Saving rate: {:.2f} images per second\".format(\n", + " float(len(exif_list)) / float((end - start).total_seconds())\n", + " )\n", + " )" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "imageprocessing", + "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.12.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/batch_processing_script.py b/batch_processing_script.py index b43cd22..fd997a9 100644 --- a/batch_processing_script.py +++ b/batch_processing_script.py @@ -3,7 +3,6 @@ import time from pathlib import Path -import numpy as np import pandas as pd from mapboxgl.utils import df_to_geojson from micasense.capture import Capture @@ -113,9 +112,7 @@ if Path("./" + warp_matrices_filename).is_file(): print("Found existing warp matrices for camera", cam_serial) - loaded = load_warp_matrices( - warp_matrices_filename, as_projective=panchro_cam - ) + loaded = load_warp_matrices(warp_matrices_filename, as_projective=panchro_cam) if panchro_cam: warp_matrices_SIFT = loaded else: diff --git a/micasense/imageutils.py b/micasense/imageutils.py index f1607a4..6229856 100644 --- a/micasense/imageutils.py +++ b/micasense/imageutils.py @@ -23,7 +23,6 @@ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import multiprocessing import os import logging diff --git a/micasense/warp_io.py b/micasense/warp_io.py index c40d77c..a3ce508 100644 --- a/micasense/warp_io.py +++ b/micasense/warp_io.py @@ -20,9 +20,7 @@ def warp_matrices_to_arrays(matrices: List[MatrixLike]) -> np.ndarray: return np.array([_matrix_to_array(m) for m in matrices], dtype=object) -def arrays_to_warp_matrices( - arrays: np.ndarray, *, as_projective: bool = True -) -> list: +def arrays_to_warp_matrices(arrays: np.ndarray, *, as_projective: bool = True) -> list: """ Convert a loaded object array back to warp matrices. @@ -44,9 +42,7 @@ def save_warp_matrices(path: Union[Path, str], matrices: List[MatrixLike]) -> No np.save(Path(path), warp_matrices_to_arrays(matrices), allow_pickle=True) -def load_warp_matrices( - path: Union[Path, str], *, as_projective: bool = True -) -> list: +def load_warp_matrices(path: Union[Path, str], *, as_projective: bool = True) -> list: """Load warp matrices written by :func:`save_warp_matrices`.""" arrays = np.load(Path(path), allow_pickle=True) return arrays_to_warp_matrices(arrays, as_projective=as_projective) diff --git a/tests/test_sift_align.py b/tests/test_sift_align.py index c2e4b4b..f02f97a 100644 --- a/tests/test_sift_align.py +++ b/tests/test_sift_align.py @@ -2,9 +2,7 @@ """Tests for SIFT_align_capture and multiprocessing spawn configuration.""" import numpy as np -import pytest -from micasense.capture import Capture from micasense.mp_config import spawn_pool diff --git a/tests/test_warp_io.py b/tests/test_warp_io.py index 7511cff..9d27135 100644 --- a/tests/test_warp_io.py +++ b/tests/test_warp_io.py @@ -5,7 +5,6 @@ from skimage.transform import ProjectiveTransform from micasense.warp_io import ( - arrays_to_warp_matrices, load_warp_matrices, save_warp_matrices, warp_matrices_to_arrays, From aa08f28cf9e1f8f4053c31944c27e90f9bc6a521 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Fri, 19 Jun 2026 08:16:18 +0200 Subject: [PATCH 8/9] fix typo in dls_irradiance call --- Batch Processing v2.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Batch Processing v2.ipynb b/Batch Processing v2.ipynb index db6f951..7b37725 100644 --- a/Batch Processing v2.ipynb +++ b/Batch Processing v2.ipynb @@ -292,7 +292,7 @@ " if panchroCam:\n", " capture.radiometric_pan_sharpened_aligned_capture(\n", " warp_matrices=warp_matrices_SIFT,\n", - " irradiance_list=capture.dls_irradiace(),\n", + " irradiance_list=capture.dls_irradiance(),\n", " img_type=img_type,\n", " write_exif=write_exif_to_individual_stacks,\n", " )\n", From 7a5bf050fc7e0371ce0250309ea4d045bfc9ffb9 Mon Sep 17 00:00:00 2001 From: Florian Maurer Date: Fri, 19 Jun 2026 09:36:16 +0200 Subject: [PATCH 9/9] update links in readme --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 5974adc..236d0b8 100644 --- a/README.md +++ b/README.md @@ -42,7 +42,7 @@ analysis. ### Tutorial Articles -[Click here to view the tutorial articles](https://micasense.github.io/imageprocessing/index.html). The set of example +[Click here to view the tutorial articles](https://brainergylab.github.io/micasense_imageprocessing/index.html). The set of example notebooks and their outputs can be viewed in your browser without downloading anything or running any code. ### How do I get set up? @@ -53,7 +53,7 @@ both before running `git clone` or you may have issues with the example data fil Next, `git clone` this repository, as it has all the code and examples you'll need. Once you have git installed and the repository cloned, you are ready to start with the first tutorial. Check out -the [setup tutorial](https://micasense.github.io/imageprocessing/MicaSense%20Image%20Processing%20Setup.html) which will +the [setup tutorial](https://brainergylab.github.io/micasense_imageprocessing/MicaSense%20Image%20Processing%20Setup.html) which will walk through installing and checking the necessary tools to run the remaining tutorials. ### MicaSense Library Usage