Library for single image processing


#1

Hi

I’m trying to process single images captured with the Sequoia mounted above a plant to capture it changing over time. Pix4D does not accept single images, unfortunately. The output I want is corrected images for each band, corrected for:

  • distortion
  • vignetting
  • reflectance

Micasense has made a library for their own sensors, and I was wondering if there is something similar for the Parrot Sequoia? This is the Micasense library in Python: https://github.com/micasense/imageprocessing

Since I couldn’t find anything so far, here is how I modified the micasense code, by changing the tags it uses (e.g. using the method in another thread for the sunshine sensor), but for vignetting, distortion and reading the QR code of the Airinov calibration plate I run into problems.

So I was hoping that someone could either point me to a working python (or other language) library or help make the micasense library work with Sequoia images. Here is the modified code of metadata.py:

#!/usr/bin/env python

import exiftool
from datetime import datetime, timedelta
import pytz
import os
import base64
import struct
import math
import pandas as pd

class Metadata(object):
    ''' Container for Sequoia image metadata'''
    def __init__(self, filename, exiftoolPath=None):
        self.xmpfile = None
        if exiftoolPath is not None:
            self.exiftoolPath = exiftoolPath
        elif os.environ.get('exiftoolpath') is not None:
            self.exiftoolPath = os.path.normpath(os.environ.get('exiftoolpath'))
        else:
            self.exiftoolPath = None
        if not os.path.isfile(filename):
            raise IOError("Input path is not a file")
        with exiftool.ExifTool(self.exiftoolPath) as exift:
            self.exif = exift.get_metadata(filename)

    def get_all(self):
        ''' Get all extracted metadata items '''
        return self.exif

    def get_item(self, item, index=None):
        ''' Get metadata item by Namespace:Parameter'''
        val = None
        try:
            val = self.exif[item]
            if index is not None:
                val = val[index]
        except KeyError:
            #print ("Item "+item+" not found")
            pass
        except IndexError:
            print("Item {0} is length {1}, index {2} is outside this range.".format(
                item,
                len(self.exif[item]),
                index))

        return val

    def size(self, item):
        '''get the size (length) of a metadata item'''
        val = self.get_item(item)
        return len(val)
    
    def print_all(self):
        for item in self.get_all():
            print("{}: {}".format(item, self.get_item(item)))

    def dls_present(self):
        return self.get_item("XMP:IrradianceList") is not None
    
    def supports_radiometric_calibration(self):
        if(self.get_item('XMP:CalibrationMeasurement')) is None:
            return False
        return True

    def position(self):
        '''get the WGS-84 latitude, longitude tuple as signed decimal degrees'''
        lat = self.get_item('EXIF:GPSLatitude')
        latref = self.get_item('EXIF:GPSLatitudeRef')
        if latref=='S':
            lat *= -1.0
        lon = self.get_item('EXIF:GPSLongitude')
        lonref = self.get_item('EXIF:GPSLongitudeRef')
        if lonref=='W':
            lon *= -1.0
        alt = self.get_item('EXIF:GPSAltitude')
        return lat, lon, alt

    def utc_time(self):
        ''' Get the timezone-aware datetime of the image capture '''
        str_time = self.get_item('EXIF:DateTimeOriginal')
        utc_time = datetime.strptime(str_time, "%Y:%m:%d %H:%M:%S")
        subsec = int(self.get_item('EXIF:SubSecTime'))
        negative = 1.0
        if subsec < 0:
            negative = -1.0
            subsec *= -1.0
        subsec = float('0.{}'.format(int(subsec)))
        subsec *= negative
        ms = subsec * 1e3
        utc_time += timedelta(milliseconds = ms)
        timezone = pytz.timezone('UTC')
        utc_time = timezone.localize(utc_time)
        return utc_time

    def dls_pose(self):
        ''' get DLS pose as local earth-fixed yaw, pitch, roll in radians '''
        yaw = float(self.get_item('XMP:Yaw')) # should be XMP.DLS.Yaw, but exiftool doesn't expose it that way
        pitch = float(self.get_item('XMP:Pitch'))
        roll = float(self.get_item('XMP:Roll'))
        return yaw, pitch, roll
    
    def dls_irradiance(self):
        # Based on http://forum.developer.parrot.com/t/details-of-irradiance-list-tag-for-sunshine-sensor-in-exif-data-of-sequoia/5261/
        # (1) Calibration data
        irradiance_calibration_measurement = self.get_item('XMP:IrradianceCalibrationMeasurement')
        calibList = [int(i) for i in irradiance_calibration_measurement.split(',')]
        relGainFac = [calibList[2]/calibList[6],
                      calibList[6]/calibList[6],
                      calibList[10]/calibList[6],
                      calibList[14]/calibList[6]]
        # (2) Measurement data
        irradiance_list = self.get_item('XMP:IrradianceList')
        irradiance_list_bytes = base64.b64decode(irradiance_list)
        irradianceData = [[]]
        for irradiance_data in struct.iter_unpack("qHHHHfff", irradiance_list_bytes):
            irradianceData.append(list(irradiance_data))
            pass
        irradianceData.pop(0)
        irradianceData = pd.DataFrame(irradianceData, columns = ['timestamp','ch0Count','ch1Count','gainMode','exposure','yaw','pitch','roll'])
        irradianceData.exposure = irradianceData.exposure/1000 # [s]
        
        #timestamp = float(irradiance_data[0])
        #ch0Count  = float(irradiance_data[1])
        #ch1Count  = float(irradiance_data[2])
        #gainMode  = int(irradiance_data[3])
        #exposure  = float(irradiance_data[4]) / 1000 # [s]
        #yaw       = float(irradiance_data[5])
        #pitch     = float(irradiance_data[6])
        #roll      = float(irradiance_data[7])
		
        # (3) Calculation of irradiance
        #gainFactor = relGainFac[gainMode]
        #irradiance = ch0Count / (gainFactor * exposure)
        # Create a list to store the data
        gainFactors = []
        for row in irradianceData['gainMode']:
            gainFactors.append(relGainFac[row])
        irradianceData['gainFactor'] = gainFactors
        irradianceData['irradiance'] = irradianceData.ch0Count / (irradianceData.gainFactor * irradianceData.exposure)
        return(irradianceData.irradiance.tail(3).mean()) # returns the mean of the last 3 irradiance values
		
    def cam_irradiance(self, p): # This does not work yet & may not belong here
        # Calculates the irradiance of any pixel p
        calibVals = self.radiometric_cal()
        exposure  = self.exposure()
        iso       = self.iso()
        aperture  = self.aperture()
        irradiance = aperture ** 2 * (p - calibVals[1]) / (calibVals[0] * exposure * iso + calibVals[2])
        return(irradiance)
		
    def reflectance(self): # This does not work yet & may not belong here
        K = 1.0 # the calibration coefficient from the reflectance panel
		# Calculate viewing angle
        yaw, pitch, roll = self.dls_pose()
        #theta = math.pi / 2 # 90° for a test
        theta = math.pi / 2 - atan(tan(roll) + tan(pitch))
		# Calculate reflectance value
        reflectanceVal = K * self.cam_irradiance()/self.dls_irradiance() / math.cos(theta)
        return(reflectanceVal)
	
    def capture_id(self):
        return self.get_item('XMP:CaptureUUID')

    def flight_id(self):
        return self.get_item('XMP:FlightUUID')

    def camera_make(self): # This is the company
        return self.get_item('EXIF:Make')

    def camera_model(self): # This is the camera model
        return self.get_item('EXIF:Model')

    def firmware_version(self):
        return self.get_item('EXIF:Software')

    def band_name(self): # A name for the image band, e.g. "Green"
        return self.get_item('XMP:BandName')
    
    def band_index(self):
        return self.get_item('XMP:RigCameraIndex')

    def exposure(self):
        return self.get_item('EXIF:ExposureTime')

    def gain(self):
        return self.get_item('EXIF:ISO')/100.0
		
    def iso(self):
        return self.get_item('EXIF:ISO')

    def aperture(self):
        return self.get_item('EXIF:ApertureValue')	

    def image_size(self):
        return self.get_item('EXIF:ImageWidth'), self.get_item('EXIF:ImageHeight')
    
    def center_wavelength(self):
        return self.get_item('XMP:CentralWavelength')

    def bandwidth(self):
        return self.get_item('XMP:WavelengthFWHM')

    def radiometric_cal(self):
        return [float(i) for i in self.get_item('XMP:SensorModel').split(',')]

    def black_level(self): # returns a mean of 4 black level measurements
        black_lvl = self.get_item('EXIF:BlackLevel').split(' ')
        total = 0.0
        num = len(black_lvl)
        for pixel in black_lvl:
            total += float(pixel)
        return total/float(num)

    #def dark_pixels(self):
    #    ''' get the average of the optically covered pixel values 
    #    Note: these pixels are raw, and have not been radiometrically
    #          corrected. Use the black_level() method for all
    #          radiomentric calibrations '''
    #    dark_pixels = self.get_item('XMP:DarkRowValue')
    #    total = 0.0
    #    num = len(dark_pixels)
    #    for pixel in dark_pixels:
    #        total += float(pixel)
    #    return total/float(num)

    def bits_per_pixel(self):
        ''' get the number of bits per pixel, which defines the maximum digital number value in the image '''
        return self.get_item('EXIF:BitsPerSample')

    def vignette_center(self):
        ''' get the vignette center in X and Y image coordinates'''
        return [(i/2) for i in self.image_size()]

    def vignette_polynomial(self):
        ''' get the radial vignette polynomial in the order it's defined in the metadata'''
        # 'XMP:VignettingPolynomial2DName' is the list of powers of 'XMP:VignettingPolynomial2D'
        return [float(i) for i in self.get_item('XMP:VignettingPolynomial2D').split(',')]

    def distortion_parameters(self):
        return [float(i) for i in self.get_item('XMP:FisheyePolynomial').split(',')]
        #nelem = self.size('XMP:FisheyePolynomial')
        #return [float(self.get_item('XMP:FisheyePolynomial', i)) for i in range(nelem)]

    def principal_point(self):
        return [float(item) for item in self.get_item('XMP:PrincipalPoint').split(',')]

    def focal_plane_resolution_px_per_mm(self):
        fp_x_resolution = float(self.get_item('EXIF:FocalPlaneXResolution'))
        fp_y_resolution = float(self.get_item('EXIF:FocalPlaneYResolution'))
        return fp_x_resolution, fp_y_resolution

    def focal_length_mm(self):
        focal_length_mm = float(self.get_item('EXIF:FocalLength'))
        return focal_length_mm

#2

@keba, I am not aware if Sequoia has any such open source library. However, you might add the functionality to code you shared above.

Vignetting correction is not hard to implement. You can find the parameters in exif data and mathematical model has already been shared here.

As for fisheye distortion and alignment, you can use this code. Code is in MATLAB but you can reprogram in python.

As far as reflectance is concerned, I would recommend you to read this document along with documents on this post.

I hope this helps. And incase you succeed in developing one library, you can share the code with the rest of community.