micasense / imageprocessing

MicaSense RedEdge and Altum image processing tutorials
https://www.micasense.com
MIT License
263 stars 152 forks source link

after running batch processing the images resulted are not aligned well #197

Closed khaledalsaih closed 1 year ago

khaledalsaih commented 1 year ago

I have used my data for batch processing to align the images using rededge camera but the images resulted are so bad, I have used this link to upload sample of the data and the results

thanks for your help.

https://drive.google.com/drive/folders/1VYdZOb5P0hk3-cYwYp5jbPFtHxboP6NX?usp=share_link

the code used as follows

from ipywidgets import FloatProgress, Layout from IPython.display import display import micasense.imageset as imageset import micasense.capture as capture import os, glob import multiprocessing import numpy as np import micasense.dls as dls import math import matplotlib.pyplot as plt

panelNames = None useDLS = True

imagePath = os.path.expanduser(os.path.join('~', 'Downloads', 'RedEdgeImageSet1', '0000SET')) panelNames = sorted(glob.glob(os.path.join(imagePath, '000', 'IMG0035*.tif'))) panelCap1 = capture.Capture.from_filelist(panelNames)

outputPath = os.path.join(imagePath, '..', 'stacks') print(outputPath) thumbnailPath = os.path.join(outputPath, '..', 'thumbnails')

overwrite = False # can be set to False to continue interrupted processing generateThumbnails = True

Define DLS sensor orientation vector relative to dls pose frame

dls_orientation_vector = np.array([0, 0, -1])

compute sun orientation and sun-sensor angles

( sun_vector_ned, # Solar vector in North-East-Down coordinates sensor_vector_ned, # DLS vector in North-East-Down coordinates sun_sensor_angle, # Angle between DLS vector and sun vector solar_elevation, # Elevation of the sun above the horizon solar_azimuth, # Azimuth (heading) of the sun ) = dls.compute_sun_angle(panelCap1.location(), panelCap1.dls_pose(), panelCap1.utc_time(), dls_orientation_vector)

Since the diffuser reflects more light at shallow angles than at steep angles,

we compute a correction for this

fresnel_correction = dls.fresnel(sun_sensor_angle)

Now we can correct the raw DLS readings and compute the irradiance on level ground

dls_irradiances = [] center_wavelengths = [] for img in panelCap1.images: dir_dif_ratio = 6.0 percent_diffuse = 1.0 / dir_dif_ratio

measured Irradiance / fresnelCorrection

sensor_irradiance = img.spectral_irradiance / fresnel_correction
untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(sun_sensor_angle))
# compute irradiance on the ground using the solar altitude angle
dls_irr = untilted_direct_irr * (percent_diffuse + np.sin(solar_elevation))
dls_irradiances.append(dls_irr)
center_wavelengths.append(img.center_wavelength)

import matplotlib.pyplot as plt

plt.scatter(center_wavelengths, dls_irradiances) plt.xlabel('Wavelength (nm)') plt.ylabel('Irradiance ($W/m^2/nm$)')

plt.show();

panelCap1.plot_undistorted_reflectance(dls_irradiances) panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] # RedEdge band_index order

panel_radiance = [0.8088253338980584,0.8246040254690193,0.7485704221000492,0.542069532328326,0.6110526527269586]

panel_radiances = np.array(panelCap1.panel_radiance()) irr_from_panel = math.pi panel_radiances / panel_reflectance_by_band dls_correction = irr_from_panel / dls_irradiances panelCap1.plot_undistorted_reflectance(dls_irradiances dls_correction)

plt.scatter(panelCap1.center_wavelengths(), panelCap1.panel_reflectance()) plt.title("Panel Reflectances") plt.xlabel("Wavelength (nm)") plt.ylabel("Reflectance")

plt.show()

panelCap = dls_correction

Allow this code to align both radiance and reflectance images; bu excluding

a definition for panelNames above, radiance images will be used

For panel images, efforts will be made to automatically extract the panel information

but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance

will need to be set in the panel_reflectance_by_band variable.

Note: radiance images will not be used to properly create NDVI/NDRE images below.

if panelNames is not None: panelCap = capture.Capture.from_filelist(panelNames) else: panelCap = None

if panelCap is not None: if panelCap.panel_albedo() is not None and not any(v is None for v in panelCap.panel_albedo()): panel_reflectance_by_band = panelCap.panel_albedo() else: panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] # RedEdge band_index order

panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
img_type = "reflectance"

else: if useDLS: img_type = 'reflectance' else: img_type = "radiance"

This progress widget is used for display of the long-running process

f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Loading") display(f) def update_f(val): if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x f.value=val

imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f) update_f(1.0)

from mapboxgl.viz import * from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between from mapboxgl.utils import create_color_stops import pandas as pd

data, columns = imgset.as_nested_lists() df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)

Insert your mapbox token here

token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw' color_property = 'dls-yaw' num_color_classes = 8

min_val = df[color_property].min() max_val = df[color_property].max()

import jenkspy breaks = jenkspy.jenks_breaks(df[color_property], n_classes=num_color_classes)

color_stops = create_color_stops(breaks,colors='YlOrRd') geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')

viz = CircleViz(geojson_data, access_token=token, color_property=color_property, color_stops=color_stops, center=[df['longitude'].median(),df['latitude'].median()], zoom=16, height='600px', style='mapbox://styles/mapbox/satellite-streets-v9') viz.show()

from numpy import array from numpy import float32

Set warp_matrices to none to align using RigRelatives

Or

Use the warp_matrices derived from the Alignment Tutorial for this RedEdge set without RigRelatives

warp_matrices = [array([[ 1.0022864e+00, -2.5218755e-03, -7.8898020e+00], [ 2.3614739e-03, 1.0036649e+00, -1.3134377e+01], [-1.7785899e-06, 1.1343118e-06, 1.0000000e+00]], dtype=float32), array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]], dtype=float32), array([[ 9.9724638e-01, -1.5535230e-03, 1.2301294e+00], [ 8.6745428e-04, 9.9738181e-01, -1.6499169e+00], [-8.2816513e-07, -3.4488804e-07, 1.0000000e+00]], dtype=float32), array([[ 1.0007139e+00, -8.4427800e-03, 1.6312805e+01], [ 6.2834378e-03, 9.9977130e-01, -1.6011697e+00], [-1.9520389e-06, -6.3762940e-07, 1.0000000e+00]], dtype=float32), array([[ 9.9284178e-01, 9.2155562e-04, 1.6069822e+01], [-3.2895457e-03, 9.9262553e-01, -5.0333548e-01], [-1.5845577e-06, -1.7680986e-06, 1.0000000e+00]], dtype=float32)]

import exiftool import datetime

This progress widget is used for display of the long-running process

f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Saving") display(f2)

def update_f2(val): f2.value = val

if not os.path.exists(outputPath): os.makedirs(outputPath) if generateThumbnails and not os.path.exists(thumbnailPath): os.makedirs(thumbnailPath)

Save out geojson data so we can open the image capture locations in our GIS

with open(os.path.join(outputPath, 'imageSet.json'), 'w') as f: f.write(str(geojson_data))

try: irradiance = panel_irradiance + [0] except NameError: irradiance = None

start = datetime.datetime.now() for i, capture in enumerate(imgset.captures): outputFilename = capture.uuid + '.tif' thumbnailFilename = capture.uuid + '.jpg' fullOutputPath = os.path.join(outputPath, outputFilename) fullThumbnailPath = os.path.join(thumbnailPath, thumbnailFilename) if (not os.path.exists(fullOutputPath)) or overwrite: if (len(capture.images) == len(imgset.captures[0].images)): capture.create_aligned_capture(irradiance_list=irradiance, warp_matrices=warp_matrices) capture.save_capture_as_stack(fullOutputPath) if generateThumbnails: capture.save_capture_as_rgb(fullThumbnailPath) capture.clear_image_data() update_f2(float(i) / float(len(imgset.captures))) update_f2(1.0) end = datetime.datetime.now()

print("Saving time: {}".format(end - start)) print("Alignment+Saving rate: {:.2f} images per second".format( float(len(imgset.captures)) / float((end - start).total_seconds())))

def decdeg2dms(dd): is_positive = dd >= 0 dd = abs(dd) minutes,seconds = divmod(dd*3600,60) degrees,minutes = divmod(minutes,60) degrees = degrees if is_positive else -degrees return (degrees,minutes,seconds)

header = "SourceFile,\ GPSDateStamp,GPSTimeStamp,\ GPSLatitude,GpsLatitudeRef,\ GPSLongitude,GPSLongitudeRef,\ GPSAltitude,GPSAltitudeRef,\ FocalLength,\ XResolution,YResolution,ResolutionUnits\n"

lines = [header] for capture in imgset.captures:

get lat,lon,alt,time

outputFilename = capture.uuid+'.tif'
fullOutputPath = os.path.join(outputPath, outputFilename)
lat,lon,alt = capture.location()
#write to csv in format:
# IMG_0199_1.tif,"33 deg 32' 9.73"" N","111 deg 51' 1.41"" W",526 m Above Sea Level
latdeg, latmin, latsec = decdeg2dms(lat)
londeg, lonmin, lonsec = decdeg2dms(lon)
latdir = 'North'
if latdeg < 0:
    latdeg = -latdeg
    latdir = 'South'
londir = 'East'
if londeg < 0:
    londeg = -londeg
    londir = 'West'
resolution = capture.images[0].focal_plane_resolution_px_per_mm

linestr = '"{}",'.format(fullOutputPath)
linestr += capture.utc_time().strftime("%Y:%m:%d,%H:%M:%S,")
linestr += '"{:d} deg {:d}\' {:.2f}"" {}",{},'.format(int(latdeg),int(latmin),latsec,latdir[0],latdir)
linestr += '"{:d} deg {:d}\' {:.2f}"" {}",{},{:.1f} m Above Sea Level,Above Sea Level,'.format(int(londeg),int(lonmin),lonsec,londir[0],londir,alt)
linestr += '{}'.format(capture.images[0].focal_length)
linestr += '{},{},mm'.format(resolution,resolution)
linestr += '\n' # when writing in text mode, the write command will convert to os.linesep
lines.append(linestr)

fullCsvPath = os.path.join(outputPath,'log.csv') with open(fullCsvPath, 'w') as csvfile: #create CSV csvfile.writelines(lines)

''' import subprocess

if os.environ.get('exiftoolpath') is not None: exiftool_cmd = os.path.normpath(os.environ.get('exiftoolpath')) else: exiftool_cmd = 'exiftool'

print(exiftool_cmd) print(fullCsvPath) print(outputPath) cmd = '{} -csv="{}" -overwrite_original="{}"'.format(exiftool_cmd, fullCsvPath, outputPath) print(cmd) subprocess.check_call(cmd) '''

khaledalsaih commented 1 year ago

closed the case cause i solved the problem