Add 2d polynomial vignertte to tutorial1

micasense · Feb 19, 2024 · b7b1649 · b7b1649
1 parent fd07d24
commit b7b1649
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 30 deletions.
diff --git 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
@@ -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