From b7b1649c0026ba3e3868a3cd2e7f206c7b399fd8 Mon Sep 17 00:00:00 2001 From: Sebastien Courroux Date: Mon, 19 Feb 2024 16:02:12 +0100 Subject: [PATCH] Add 2d polynomial vignertte to tutorial1 --- MicaSense Image Processing Tutorial 1.ipynb | 24 ++++++--- micasense/utils.py | 58 +++++++++++++-------- 2 files changed, 52 insertions(+), 30 deletions(-) diff --git a/MicaSense Image Processing Tutorial 1.ipynb b/MicaSense Image Processing Tutorial 1.ipynb index 3b9ab686..5d5aef25 100644 --- a/MicaSense Image Processing Tutorial 1.ipynb +++ b/MicaSense Image Processing Tutorial 1.ipynb @@ -29,11 +29,24 @@ "%matplotlib inline\n", "\n", "imagePath = os.path.join('.','data','REDEDGE-MX')\n", - "imageName = os.path.join(imagePath,'IMG_0001_4.tif')\n", + "panelImageName = os.path.join(imagePath,'IMG_0001_4.tif')\n", + "ulx = 610 # upper left column (x coordinate) of panel area\n", + "uly = 610 # upper left row (y coordinate) of panel area\n", + "lrx = 770 # lower right column (x coordinate) of panel area\n", + "lry = 760 # lower right row (y coordinate) of panel area\n", + "flightImageName = os.path.join(imagePath,'IMG_0020_4.tif')\n", + "\n", + "# imagePath = os.path.join('.','data','REDEDGE-P')\n", + "# panelImageName = os.path.join(imagePath,'IMG_0000_4.tif')\n", + "# ulx = 520 # upper left column (x coordinate) of panel area\n", + "# uly = 505 # upper left row (y coordinate) of panel area\n", + "# lrx = 670 # lower right column (x coordinate) of panel area\n", + "# lry = 660 # lower right row (y coordinate) of panel area\n", + "# flightImageName = os.path.join(imagePath,'IMG_0011_4.tif')\n", "\n", "# Read raw image DN values\n", "# reads 16 bit tif - converts 12 bit to 16 bit if necessary \n", - "imageRaw=plt.imread(imageName)\n", + "imageRaw=plt.imread(panelImageName)\n", "\n", "# Display the image\n", "fig, ax = plt.subplots(figsize=(8,6))\n", @@ -91,7 +104,7 @@ "if os.name == 'nt':\n", " exiftoolPath = os.environ.get('exiftoolpath')\n", "# get image metadata\n", - "meta = metadata.Metadata(imageName, exiftool_path=exiftoolPath)\n", + "meta = metadata.Metadata(panelImageName, exiftool_path=exiftoolPath)\n", "cameraMake = meta.get_item('EXIF:Make')\n", "cameraModel = meta.get_item('EXIF:Model')\n", "firmwareVersion = meta.get_item('EXIF:Software')\n", @@ -200,10 +213,6 @@ "outputs": [], "source": [ "markedImg = radianceImage.copy()\n", - "ulx = 610 # upper left column (x coordinate) of panel area\n", - "uly = 610 # upper left row (y coordinate) of panel area\n", - "lrx = 770 # lower right column (x coordinate) of panel area\n", - "lry = 760 # lower right row (y coordinate) of panel area\n", "cv2.rectangle(markedImg,(ulx,uly),(lrx,lry),(0,255,0),3)\n", "\n", "# Our panel calibration by band (from MicaSense for our specific panel)\n", @@ -293,7 +302,6 @@ "metadata": {}, "outputs": [], "source": [ - "flightImageName = os.path.join(imagePath,'IMG_0001_4.tif')\n", "flightImageRaw=plt.imread(flightImageName)\n", "plotutils.plotwithcolorbar(flightImageRaw, 'Raw Image')\n", "\n", diff --git a/micasense/utils.py b/micasense/utils.py index 65fec77d..a6f83981 100755 --- a/micasense/utils.py +++ b/micasense/utils.py @@ -77,33 +77,47 @@ def raw_image_to_radiance(meta, image_raw): def vignette_map(meta, x_dim, y_dim): - # get vignette center - x_vignette = float(meta.get_item('XMP:VignettingCenter', 0)) - y_vignette = float(meta.get_item('XMP:VignettingCenter', 1)) - - # get vignette polynomial - nvignette_poly = meta.size('XMP:VignettingPolynomial') - vignette_poly_list = [float(meta.get_item('XMP:VignettingPolynomial', i)) for i in range(nvignette_poly)] - - # reverse list and append 1., so that we can call with numpy polyval - vignette_poly_list.reverse() - vignette_poly_list.append(1.) - vignette_poly = np.array(vignette_poly_list) - - # perform vignette correction - # get coordinate grid across image + # get coordinate grid across image, seem swapped because of transposed vignette x, y = np.meshgrid(np.arange(x_dim), np.arange(y_dim)) - # meshgrid returns transposed arrays x = x.T y = y.T - # compute matrix of distances from image center - r = np.hypot((x - x_vignette), (y - y_vignette)) - - # compute the vignette polynomial for each distance - we divide by the polynomial so that the - # corrected image is image_corrected = image_original * vignetteCorrection - vignette = 1. / np.polyval(vignette_poly, r) + vignetting_center_size = meta.size('XMP:VignettingCenter') + vignetting_polynomial_2d_size = meta.size('XMP:VignettingPolynomial2D') + + # if we have a radial poly + if vignetting_center_size > 0: + # get vignette center + x_vignette = float(meta.get_item('XMP:VignettingCenter', 0)) + y_vignette = float(meta.get_item('XMP:VignettingCenter', 1)) + + # get vignette polynomial + nvignette_poly = meta.size('XMP:VignettingPolynomial') + vignette_poly_list = [float(meta.get_item('XMP:VignettingPolynomial', i)) for i in range(nvignette_poly)] + + # reverse list and append 1., so that we can call with numpy polyval + vignette_poly_list.reverse() + vignette_poly_list.append(1.) + vignette_poly = np.array(vignette_poly_list) + + # compute matrix of distances from image center + r = np.hypot((x - x_vignette), (y - y_vignette)) + + # compute the vignette polynomial for each distance - we divide by the polynomial so that the + # corrected image is image_corrected = image_original * vignetteCorrection + vignette = 1. / np.polyval(vignette_poly, r) + elif vignetting_polynomial_2d_size > 0: # 2d polynomial + xv = x.T / x_dim + yv = y.T / y_dim + k = meta.vignette_polynomial2D() + e = meta.vignette_polynomial2Dexponents() + p2 = np.zeros_like(xv, dtype=float) + for i, c in enumerate(k): + ex = e[2 * i] + ey = e[2 * i + 1] + p2 += c * xv ** ex * yv ** ey + vignette = (1. / p2).T return vignette, x, y