Zombie Door

Run!!!

Zombies are pretty cool.  This post describes something a little less cool, but uses zombies to explain the concept (in a shallow, transparent attempt to capture your attention!)

Zombie Door Method

Imagine we want to generate a uniformly distributed random sampling in some complex space that our random number generator does not directly support.

Let me start with an simple example.  Imagine we have a random number generator that produces a random integer between 1 and 100.  However, we actually want to generate random numbers between 41 and 50.  (I know there are better ways to do this, but stick with me for this example.)  Let’s solve this with the zombie door method.

  • Build a wall and label it 1-100 where 1 is to the far left, and 100 is to the far right.
  • Cut a door in your wall from 41 to 50.
  • Now create random zombies starting anywhere from 1 to 100 and let them walk straight towards the wall.
  • The zombies will lurch across the room and eventually hit wall, explode, and dissolve in a fizzing puddle of bloody goo … or whatever it is that zombies do when they die.
  • The zombies that are lucky enough to walk through the door survive!

For this example it would be easy to scale and offset the usually available random number generator that produces a floating point value between 0 and 1, but it illustrates the basic approach of the ‘zombie door’ method.

A More Interesting Example

Imagine we have an arbitrary polygon outline and we want to splatter it with random circles.  However, we only want circles that are completely within the polygon.  (And we want our random circles to be a true  unbiased, uniformly distributed random sample.)  This example is just like the simple ‘wall’ example except now we have gone ‘2D’.

Imagine an arbitrary polygon shape:

We would like to fill this shape with 200 random circles, making sure none of our circles straddle the boundary or lie outside the polygon.  We want an even distribution over the interior area of our polygon.

We can do this by generating random circles within the min/max range of our shape, and then testing if they lie completely inside the polygon.  We reject any outliers and keep the inliers.  It’s very simple, but very cool because … you know … zombies!  Here is a picture of how our circles turned out:

If you are at all curious as to what became of our dead zombies (all the ones that splatted against the wall) here is the picture of that:

If you the reader are concerned that I am making light of a very real pending zombie apocalypse issue then here is my best tip for protecting your residence:

Finally, here is the python script I used to generate the pictures in this posting.  I leverage the very, very, very cool GPC polygon clipping library to test if the outline polygon ‘covers’ the circle.

#!/usr/bin/python
 
import random
 
# polygon (GPC) python package: https://pypi.python.org/pypi/Polygon2
import Polygon
import Polygon.IO
import Polygon.Shapes
 
# this will seed the random number generator with current time and
# force different results on each run.
random.seed()
 
# this is the bounds of the shape that will contain our splatters
outline = [ [0.5, -0.01], [1.01, -0.01], [0.5, 1.01], [-0.01, 1.01] ]
 
# define circle properties
num_circles = 200
min_r = 0.002
max_r = 0.02
 
# no closer than this to the boundary
margin = 0
#margin = 0.002
 
# create the polygon template (outline)
template = Polygon.Polygon(outline)
 
# make the shape a bit more interesting
c = Polygon.Shapes.Circle(radius=0.4, center=(1, 0.5), points=32)
template = template - c
c = Polygon.Shapes.Circle(radius=0.3, center=(0, 1), points=32)
template = template - c
 
# determine max/min of template
min_x = max_x = outline[0][0]
min_y = max_y = outline[0][1]
for p in outline:
 if p[0] < min_x: min_x = p[0]
 if p[0] > max_x: max_x = p[0]
 if p[1] < min_y: min_y = p[1]
 if p[1] > max_y: max_y = p[1]
 
print 'template bounds:', min_x, min_y, 'to', max_x, max_y
print 'radius range:', min_r, max_r
print 'margin:', margin
print 'num circles:', num_circles
 
# generate splats using zombie door method
circles = []
discarded = []
while len(circles) < num_circles:
 x = random.uniform(min_x, max_x)
 y = random.uniform(min_y, max_y)
 r = random.uniform(min_r, max_r)
 
 # make the circle
 c = Polygon.Shapes.Circle(radius=r, center=(x, y), points=32)
 
 # make the circle padded with extra margin
 cm = Polygon.Shapes.Circle(radius=(r+margin), center=(x, y), points=32)
 
 if template.covers(cm):
  # circle + margin fully contained inside the template
  circles.append(c)
 else:
  discarded.append(c)
 
# assemble final polygons and write output
final = Polygon.Polygon()
for c in circles:
 final += c
Polygon.IO.writeGnuplot('in.plt', [template, final])
Polygon.IO.writeSVG('in.svg', [final], fill_color=(0,0,0))
 
reject = Polygon.Polygon()
for c in discarded:
 reject += c
Polygon.IO.writeGnuplot('out.plt', [template, reject])
Polygon.IO.writeSVG('out.svg', [reject], fill_color=(0,0,0))

Failure is not fatal

This post is penned during a moment of extreme frustration, beware!

Kobayashi Maru

https://en.wikipedia.org/wiki/Kobayashi_Maru

One of the reasons I loved the original Star Trek series is because no matter what the odds, no matter how hopeless the circumstances, no matter how impossible the foe, Captain Kirk always found a way to think his way out of the mess.  He never ultimately failed or lost to an opponent, not once, not ever.  That makes a great hero and fun TV!  Fictional super heroes do things that normal human beings could never possibly do … like fly, or be stronger than steel, or always win.

Stress and the Brain

I don’t have time to read stuff like the following link, especially when I’m coming up short of a promised deadline.  Maybe you do?  http://www.health.harvard.edu/staying-healthy/understanding-the-stress-response

I’m told that when we begin to get stressed, the front area of the brain that is responsible for logic and reason starts to shut down, and command functions begin to be transferred back to the “fight or flight” portion of the brain.  I think about standing up in front of a group and speaking, then sitting down and wondering what I even said?  I think about arguments that got out of hand?  Where was the front part of my brain in all of that?  I think about looming deadlines and mounting stress … and … and … and mounting stress!

Recursive Stress

My job largely amounts to puzzle solving.   I love the process and I love finding clever solutions.  But if I ask you a riddle or give you a logic problem, can you give me a specific estimate of how much time it will take you to solve it?  That’s not how puzzle solving works, it’s not a step by step recipe that leads to a solution in a known time.  Failing to solve the problem in time stresses me out!  What is needed in these situations is clear, logical, and calm thinking.  But that is the first part of the brain to turn off during stressful situations!  It’s exactly the part of the brain we desperately need the most.  I know all this, and I watch helplessly as it happens.  What does that create?  More stress of course which accelerates the process of losing the most important part of my brain!

What is the solution?

No, seriously, what is the solution???

People often say they do their best work under pressure.  I know for myself, I do my worst work under pressure.  I strive whenever possible to get a long head start on a complex and difficult task.  I strive whenever possible to identify and solve the hardest parts of the task first.  But that isn’t always possible.

So instead I sometimes see failure coming weeks away, maybe like an asteroid on a collision course with earth.  I’m very serious about the task, I do everything I possibly can, I pour in all my energy and expertise, but it’s not always enough.  Things I thought would be easy turn out to be 10x more difficult than imagined.  Things that were working break for unexpected reasons.  Things that shouldn’t take time, take way too much precious time.

Captain Kirk to the rescue?  Sadly, no … he is a fictional character.  In the real world the asteroid looms bigger and bigger, it’s trajectory is a mathematical certainty.  The physics of the impact can largely be predicted.  At some point it becomes clear my efforts will fall short and there’s nothing left to do but watch the asteroid in it’s last few hours of flight.  Then <boom>.

Is it just me that fails colossally?

It usually seems like I’m the only one that makes a miserable mess of things I try to do, the things I’ve promised I could do, things I’ve been paid to do.  Everyone else is posting about their giant success on facebook.  Everyone else’s resume is a spotless collection of triumphs.  But not me.  Maybe once or twice I got lucky and the rest of the time is a huge struggle!  Honestly though, the only reason I’m posting this is because I know it’s not just me.  Any sports fans out there?  How many teams and players succeed to win the championship every season?  What percentage of players ever win a championship in their whole career?  Political success and failure?  How many new businesses succeed versus failing?

High Profile Failures

By mentioning specific companies, I don’t mean to imply specific people or imply anything negative here.  My intent here is to show we are all in this together and we all, even the best and most successful of us, suffer set backs in our work.  I live and work in the world of drones (or small unmanned aerial systems, aka UAS’s).  This is a tough business.  For all the hype and excitement even big companies can struggle.  Gopro recently did a full recall of their new drone product.  Hopefully they’ll try again in 2017, and hopefully the process will go better for them.  Recently 3DR, the king of DIY drones announced they were cancelling all their hardware efforts to focus on a software package for managing data collected by drones.  Parrot (another big name in the small drone market) just announced layoffs.  Edit: 12 Jan, Lily just announces it is dropping out of the drone race and shutting down.  Edit: Facebook Aquila crashed on first test flight.  Edit:  Titan (google’s own high altitude effort solar powered effort) is shut down.  It’s tough, even for the big guys with enough money to hire the best engineers, best managers and do professional marketing.

There are even higher profile failures than these … Chernobyl, TIger Woods, the Titanic, Exon Valdez, the Challenger, and most Sylvester Stallone movies except for Rocky I.

The Aftermath

So the asteroid hit.  In the last moments we just threw up our hands, gave up, and just watched it come in and do it’s destruction.  The dust is settling, what happens next?  Maybe the asteroid wasn’t as big as we imagined?  Maybe the damage not as severe?  Maybe life goes on.  In a work context, maybe you take a big professional hit on your reputation?  Maybe you don’t?  Maybe it’s about survival and living to fight another day?

Failures suck.  Failures happen to everyone.  Failures are [usually] not fatal.  The sun will still rise tomorrow for most of us.

Survival – Yes?!?

If you are reading this, you are still here and still surviving.  That’s great news!  Hopefully I’m here too!  Lets all live to fight another day.  Let’s all help each other out when we can!  There is a word called “grace” which is worth looking up if you don’t know what it means.  It’s a quantity that we all need more of and we all need to give each other in big healthy doses!

“Success is not final; failure is not fatal. It is the courage to continue that counts.” — Budweiser ad, 1938.  (Not Winston Churchill)

Synthetic Air Data (an afternoon hack)

Motivation

Air France #447
Air France #447

On June 1, 2009 Air France flight #447 disappeared over the Atlantic Ocean.  The subsequent investigation concluded “that the aircraft crashed after temporary inconsistencies between the airspeed measurements – likely due to the aircraft’s pitot tubes being obstructed by ice crystals – caused the autopilot to disconnect, after which the crew reacted incorrectly and ultimately caused the aircraft to enter an aerodynamic stall from which it did not recover.”  https://en.wikipedia.org/wiki/Air_France_Flight_447

This incident along with a wide variety of in-flight pitot tube problems across the aviation world have led the industry to be interested in so called “synthetic airspeed” sensors.  In other words, is it possible to estimate the aircraft’s airspeed by using a combination of other sensors and techniques when we are unable to directly measure airspeed with a pitot tube?

Basic Aerodynamics

Next, I need to say a quick word about basic aerodynamics with much hand waving.  If we fix the elevator position of an aircraft, fix the throttle position, and hold a level (or fixed) bank angle, the aircraft will typically respond with a slowly damped phugoid and eventually settle out at some ‘trimmed’ airspeed.

If the current airspeed is faster than the trimmed airspeed, the aircraft will have a positive pitch up rate which will lead to a reduction in airspeed.  If the airspeed is slower than the trimmed airspeed, the aircraft will have a negative pitch rate which will lead to an acceleration.  https://en.wikipedia.org/wiki/Phugoid

The important point is that these variables are somehow all interrelated.  If you hold everything else fixed, there is a distinct relationship between airspeed and pitch rate, but this relationship is highly dependent on the current position of elevator, possibly throttle, and bank angle.

Measurement Variables and Sensors

In a small UAS under normal operating conditions, we can measure a variety of variables with fairly good accuracy.  The variables that I wish to consider for this synthetic airspeed experiment are: bank angle, throttle position, elevator position, pitch rate, and indicated airspeed.

We can conduct a flight and record a time history of all these variables.  We presume that they have some fixed relationship based on the physics and flight qualities of the specific aircraft in it’s current configuration.

It would be possible to imagine some well crafted physics based equation that expressed the true relationship between these variables …. but this is a quick afternoon hack and that would require too much time and too much thinking!

Radial Basis Functions

Enter radial basis functions.  You can read all about them here:  https://en.wikipedia.org/wiki/Radial_basis_function

From a practical perspective, I don’t really need to understand how radial basis functions work.  I can simply write a python script that imports the scipy.interpolate.Rbf module and just use it like a black box.  After that, I might be tempted to write a blog post, reference radial basis functions, link to wikipedia, and try to sound really smart!

Training the Interpolater

Step one is to dump the time history of these 5 selected variables into the Rbf module so it can do it’s magic.  There is a slight catch, however.  Internally the rbf module creates an x by x matrix where x is the number of samples you provide.   With just a few minutes of data you can quickly blow up all the memory on your PC.  As a work around I split the entire range of all the variables into bins of size n.  In this case I have 4 independent variables (bank angle, throttle position, elevator position, and pitch rate) which leads to an nnnn matrix.  For dimensions in the range of 10-25 this is quite manageable.

Each element of the 4 dimensional matrix becomes a bin that holds he average airspeed for all the measurements that fall within that bin.  This matrix is sparse, so I can extract just the non-zero bins (where we have measurement data) and pass that to the Rbf module.  This accomplishes two nice results: (1) reduces the memory requirements to something that is manageable, and (2) averages out the individual noisy airspeed measurements.

Testing the Interpolater

Now comes the big moment!  In flight we can still sense bank angle, throttle position, elevator position, and pitch rate.  Can we feed these into the Rbf interpolater and get back out an accurate estimate of the airspeed?

Here is an example of one flight that shows this technique actually can produce some sensible results.  Would this be close enough (with some smoothing) to safely navigate an aircraft through the sky in the event of a pitot tube failure?  Could this be used to detect pitot tube failures?  Would this allow the pitot tube to be completely removed (after the interpolater is trained of course)?

Synthetic airspeed estimate versus measured airspeed.
Zoom on one portion of the flight.
Zoom in on another portion of the flight.

Source Code

The source code for this experimental afternoon hack can be found here (along with quite a bit of companion code to estimate aircraft attitude and winds via a variety of algorithms.)

https://github.com/AuraUAS/navigation/blob/master/scripts/synth_asi.py

Conclusions

This is the results of a quick afternoon experiment.  Hopefully I have showed that creating a useful synthetic airspeed sensor is possible.  There are many other (probably better) ways a synthetic air speed sensor could be derived and implemented.  Are there other important flight variables that should be considered?  How would you create an equation that models the physical relationship between these sensor variables?  What are your thoughts?

Howto: Action Cam HUD overlay

NOTICE: This is a draft document and in the process of being written.  It is incomplete and subject to change at any time without notice.  

In this post I share my process and tools for creating HUD overlays on a flight video.  Here is a quick overview of the process: (1) Calibrate your camera, (2) Extract the roll rate information from the flight video, (3) Automatically and perfectly time correlate the video with the flight data, (4) Render a new video with the HUD overlay.

Source Code

Please be aware that this code has not yet gone through a v1.0 release process or outside review, so you are likely to run into missing python packages or other unanticipated issues.  I hope that a few brave souls will plow through this for themselves and help me resolve poor documentation or gaps in the package requirements.  All the code is open-source (MIT license) and available here:

https://github.com/UASLab/ImageAnalysis

Camera Calibration

Let’s jump right into it.  The very first thing that needs to be done is calibrate your specific camera.  This might sound difficult and mysterious, but have no fear, there is a script for everything!  Every camera (especially cheaper action cameras) has unique lens imperfections.  It is best to do your own calibration for each of your cameras rather than trying save some time and copy someone else’s configuration.

The camera calibration process involves feeding several images of a checkerboard calibration pattern to a calibration script.  Each image is analyzed to locate the checkerboard pattern.  This information is then passed to a solver that will find your camera’s specific calibration matrix and lens distortion parameters.

To make this process easier, the calibration script expects a short 1-2 minute movie.  It processes each frame of the movie, locates the checkboard pattern and stashes that information away.  After all frames are processed, the script samples ‘n’ frames from the movie (where ‘n’ is a number around 100-200) and uses those frames to solve for your camera’s calibration.  The reason that all frames are not used is because when ‘n’ starts pushing 250-300+, the solver begins to take a long time where long is measured in hours not minutes.

Here is a sample calibration video.  The goal is to always keep the checkerboard pattern fully in view while moving it around, closer, further, to different parts of the frame, and from different offset angles.  It is also important to hold the camera steady and move it slowly to avoid the effects of blurring and rolling shutter.

Now save your movie to your computer and run:

calibrate_movie.py --movie <name>

The script will run for quite some time (be patient!) but will eventually spit out a camera calibration matrix and set of lens distortion parameters. Save these somewhere (copy/paste is your friend.)

Extract the Roll Rate Information from the Flight Video

First take the camera calibration and distortion parameters derived in step one, and copy them into the gyro rate estimation script.

Next, run the script.  This script will detect features in each frame, match the found features with the previous frame, and compute the amount of rotation and translation from each frame to the next.  (It also is a hacky image stabilization toy).

1-est-gyro-rates.py --scale 0.4 --movie <name> --no-equalize

Here is an example of the script running and detecting features.  Notice that a very significant portion of the frame is covered by the aircraft nose and the prop covers most of the remaining area.  That’s ok!  The gyro estimation is still good enough to find the correlation with the flight data.

Correlate the Video with the Flight Data

The next script actually performs both of the last two steps (correlation and rendering.)

2-gen-hud-overlay.py --movie <name> --aura-dir <flight data dir> --resample-hz 30 --scale 0.45

The script loads the flight data log and the movie data log (created by the previous script).  It resamples them both at a common fixed rate (30hz in the above example.)  Then it computes the best correlation (or time offset) between the two.

Here is a plot of the roll rate estimate from the video overlaid with the actual roll gyro.  You can see there are many differences, but the overall trend and pattern leads to only one possible correct time correlation.

video-flight-log-correlation

This graph is especially noisy because only a large portion of the outside view is obscured by the nose and the visible portion is further obscured by the propeller.  But it is ok, the correlation process is magic and is really good at finding the best true correlation.  The next plot shows the results we can attain when we have more idea conditions with an unobstructed view.  Here is a plot that shows the video roll estimate and the actual roll gyro are almost perfectly in agreement.

overlay

Taking a step back, what did we just do there?  Essentially, we have created an automated way to align the video frames with the flight data log.  In other words, for any video frame number, I can compute the exact time in the flight log, and for any point in the flight log, I can compute the corresponding video frame number.  Now all that is left is to draw the exact current flight data (that we now have a way to find) on top of the video.

I look forward to your comments and questions!

This tutorial is far from complete and I know there are some built in assumptions about my own system and aircraft cooked into the scripts.  Please let me know your questions or experiences and I will do my best to answer or improve the code as needed.

Autopilot Visualization

Blending real video with synthetic data yields a powerful and cool! way to visualize your kalman filter (attitude estimate) as well as your autopilot flight controller.

screenshot-from-2016-11-30-14-31-14

Conformal HUD Elements

Conformal definition: of, relating to, or noting a map or transformation in which angles and scale are preserved.  For a HUD, this means the synthetic element is drawn in a way that visually aligns with the real world.  For example: the horizon line is conformal if it aligns with the real horizon line in the video.

  • Horizon line annotated with compass points.
  • Pitch ladder.
  • Location of nearby airports.
  • Location of sun, moon, and own shadow.
  • If alpha/beta data is avaliable, a flight path marker is drawn.
  • Aircraft nose (i.e. exactly where the aircraft is pointing towards.)

Nonconformal HUD Elements

  • Speed tape.
  • Altitude tape.
  • Pilot or autopilot ‘stick’ commands.

Autopilot HUD Elements

  • Flight director vbars (magenta).  These show the target roll and pitch angles commanded by the autopilot.
  • Bird (yellow).  This shows the actual roll and pitch of the aircraft.  The autopilot attempts to keep the bird aligned with the flight director using aileron and elevator commands.
  • Target ground course bug (show on the horizon line) and actual ground course.
  • Target airspeed (drawn on the speed tape.)
  • Target altitude (drawn on the altitude tape.)
  • Flight time (for referencing the flight data.)

Case Study #1: EKF Visualization

(Note this video was produced earlier in the development process and doesn’t contain all the HUD elements described above.)

What to watch for:

  • Notice the jumpiness of the yellow “v” on the horizon line.  This “v” shows the current estimated ground track, but the jumpiness points to an EKF tuning parameter issue that has since been resolved.
  • Notice a full autonomous wheeled take off at the beginning of the video.
  • Notice some jumpiness in the HUD horizon and attitude and heading of the aircraft.  This again relates back to an EKF tuning issue.

I may never have noticed the EKF tuning problems had it not been for this visualization tool.

Case Study #2: Spin Testing

What to watch for:

  • Notice the flight path marker that shows actual alpha/beta as recorded by actual alpha/beta airdata vanes.
  • Notice how the conformal alignment of the hud diverges from the real horizon especially during aggressive turns and spins.  The EKF fits the aircraft attitude estimate through gps position and velocity and aggressive maneuvers lead to gps errors (satellites go in and out of visibility, etc.)
  • Notice that no autopilot symbology is drawn because the entire flight is flown manually.

Case Study #3: Skywalker Autopilot

What to watch for:

  • Notice the yellow “v” on the horizon is still very jumpy.  This is the horizontal velocity vector direction which is noisy due to EKF tuning issues that were not identified and resolved when this video was created.  In fact it was this flight where the issue was first noticed.
  • Notice the magenta flight director is overly jumpy in response to the horizontal velocity vector being jumpy.  Every jump changes the current heading error which leads to a change in roll command which the autopilot then has to chase.
  • Notice the flight attitude is much smoother than the above Senior Telemaster flight.  This is because the skywalker EKF incorporates magnetometer measurements as well as gps measurements and this helps stabilize the filter even with poor noise/tuning values.
  • You may notice some crazy control overshoot on final approach.  Ignore this!  I was testing an idea and got it horribly wrong.  I’m actually surprised the landing completed successfully, but I’ll take it.
  • Notice in this video the horizon stays attached pretty well.  Much better than in the spin-testing video due to the non-aggressive flight maneuvers, and much better than the telemaster video due to using a more accurate gps: ublox7p versus ublox6.  Going forward I will be moving to the ublox8.

Case Study #4: Results of Tuning

What to watch for:

  • Notice that visually, the HUD horizon lines stays pegged to the camera horizon within about a degree for most of this video.  The EKF math says +/-3 standard deviations is about 1.4 degrees in pitch and roll.
  • You may notice a little more variation in heading.  +/-3 standard deviations in heading is roughly 4 degrees.
  • Now that the EKF is tamed a bit better, we can start to tune the PID’s and go after some more subtle improvements in flight control.  For example, this skywalker is kind of a floppy piece of foam.  I estimate that I have to hold a 4-5 degree right bank to fly straight.  We can begin to account for these aircraft specific nuances to improve tracking, autoland performance, etc.

Howto?

These flight visualization videos are created with an automated process using open source tools and scripts. I have started a post on how to create these videos yourself.

Mistakes!

dj-mistakes-homer-dohmistake-pano_13891

I make thousands of mistakes a day, mistakes typing, mistakes coding software, mistakes driving, mistakes walking, forgetting to order my sandwich without mayo, etc.  Most of the time they are immediately obvious — a red squiggly line under a word I mistyped, a compiler spewing an error message on line #42, a stubbed toe, my gps suggesting a u-turn at the next intersection, etc.

mistakes

But what happens when the mistake isn’t obvious, isn’t noticed immediately, and doesn’t cause everything around me to immediately fail?  Often these mistakes can have a long lifespan.  Often we discover them when we are looking for something else.

Mistakes from the Trenches.

I wanted to write about a few subtle unnoticed mistakes that lurked in the AuraUAS code for quite some time.

Temperature Calibration #Fail

AuraUAS has a really cool capability where it can estimate the bias (error) of the accelerometers during flight.  The 15-state EKF does this as part of it’s larger task of estimating the aircraft’s attitude, location, and velocity.  These bias estimates along with the corresponding IMU temperature can be used to build up a temperature calibration fit for each specific IMU based on flight data over time.  The more you fly in different temperature conditions, the better your temperature calibration becomes.  Sweet!  Calibrated accelerometers are important because accel calibration errors directly translate to errors in initial roll and pitch estimates (like during launch or take off where these values can be critical.)  Ok, the EKF will sort them out once in the air, because that is a cool feature of the EKF, but it can’t work out the errors until after flying a bit.

The bias estimates and temperature calibration fit are handled by post-flight python scripts that work with the logged flight data.  Question: should I log the raw accel values or should I log the calibrated accel values.  I decided I should log the calibrated values and then use the inverse calibration fit function to derive the original raw values after the flight.  Then I use these raw values to estimate the bias (errors), add the new data to the total collection of data for this particular IMU, and revise the calibration fit.  The most straightforward path is to log calibrated values on board during flight (in real time) and push the complicated stuff off into post processing.

However, I made a slight typo in the property name of the temperature range limits for the fit (we only fit within the range of temperatures we have flight data for.)  This means the on-board accel correction was forcing the temperature to 27C (ignoring the actual IMU temperature.)  However, when backing out the raw values in post processing, I was using the correct IMU temperature and thus arriving at a wrong raw value.  What a mess.  That means a year of calibration flight data is basically useless and I have to start all my IMU calibration learning over from scratch.  So I fixed the problem and we go forward from here with future flights producing a correct calibration.

Integer Mapping #Fail

This one is subtle.  It didn’t produce incorrect values, it simply reduced the resolution of the IMU gyros by a factor of 4 and the accels by a factor of 2.

Years ago when I first created the apm2-sensor firmware — that converts a stock APM2 (atmega2560) board into a pure sensor head — I decide to change the configured range of the gyros and accels.  Instead of +/-2000 degrees per second, I set the gyros for +/-500 degrees per second.  Instead of +/-8 g’s on the accels, I set them for +/- 4 g’s.  The sensed values get mapped to a 16 bit integer, so using a smaller range results in more resolution.

The APM2 reads the raw 16 bit integer values from the IMU and converts this to radians per second.  However, when the APM2 sends these values to the host, it re-encodes them from a 4-byte float to a 2-byte (16-bit) integer to conserve bandwidth.  Essentially this undoes the original decoding operation to efficiently transmit the values to the host system.  The host reads the encoded integer value and reconverts it into radians per second for the gyros (or mps^2 for the accels.)

The problem was that for encoding and decoding between the APM2 and the host, I used the original scaling factor for +/-2000 dps and +/-8g, not the correct scaling factor for the new range I had configured.  This mistake caused me to lose all the resolution I intended to gain.  Because the system produced the correct values on the other end, I didn’t notice this problem until someone asked me exactly what resolution the system produced, which sent me digging under the hood to refresh my memory.

This is now fixed in apm2-sensors v2.52, but requires a change to the host software as well so the encoding and decoding math agrees.  Now the IMU reports the gyro rates with a resolution of 0.015 degrees per second where as previously the resolution was 0.061 degrees per second.  Both are actually pretty good, but it pained me to discover I was throwing away resolution needlessly.

Timing #Fail

This one is also very subtle; timing issues often are.  In the architecture of the AuraUAS flight controller there is an APM2 spitting out new sensor data at precisely 100 hz.  The host is a beaglebone (or any linux computer) running it’s own precise 100 hz main loop.  The whole system runs at 100 hz throughput and life is great — or so I thought.

I had been logging flight data at 25hz which has always been fine for my own needs.  But recently I had a request to log the flight data at the full 100 hz rate.  Could the beaglebone handle this?  The answer is yes, of course, and without any trouble at all.

A question came up about logging high rate data on the well known PX4, so we had a student configure the PX4 for different rates and then plot out the time slice for each sample.  We were surprised at the huge variations in the data intervals, ranging from way too fast, to way too slow, and rarely exactly what we asked for.

I know that the AuraUAS system runs at exactly 100hz because I’ve been very careful to design it that way.  Somewhat smugly I pulled up a 100hz data set and plotted out the time intervals for each IMU record.  The plot surprised me — my timings were all over the map and not much better than the PX4.  What was going on?

I took a closer look at the IMU records and noticed something interesting.  Even though my main loop was running precisely and consistently at 100 hz, it appeared that my system was often skipping every other IMU record.  AuraUAS is designed to read whatever sensor data is available at the start of each main loop iteration and then jump into the remaining processing steps.  Because the APM2 runs it’s own loop timing separate from the host linux system, the timing between sending and receiving (and uart transferring) can be misalligned so that when the host is ready to read sensor data, there might not be any yet, and next time there may be 2 records waiting.  It is subtle, but communication between to free running processor loops can lead to issues like this.  The end result is usually still ok, the EKF handles variable dt just fine, the average processing rate maybe drops to 50hz, and that’s still just fine for flying an airplane around the sky … no big deal right?  And it’s really not that big of a deal for getting the airplane from point A to point B, but if you want to do some analysis of the flight data and want high resolution, then you do have a big problem.

What is the fix?  There are many ways to handle timing issues in threaded and distributed systems.  But you have to be very careful, often what you get out of your system is not what you expected or intended.  In this case I have amended my host system’s main loop structure to throw away it’s own free running main loop.  I have modified the APM2 data output routine to send the IMU packet at the end of each frame’s output to mark the end of data.  Now the main loop on the host system reads sensor data until it receives an IMU packet.  Then and only then does it drop through to the remaining processing steps.  This way the timing of the system is controlled precisely by the APM2, the host system’s main loop logic is greatly simplified, and the per frame timing is far more consistent … but not consistent enough.

The second thing I did was to include the APM2 timestamp with each IMU record.  This is a very stable, consistent, and accurate timestamp, but it counts up from a different starting point than the host.  On the host side I can measure the difference between host clock and APM2 clock, low pass filter the difference and add this filtered difference back to the APM2 timestamp.  The result is a pretty consistent value in the host’s frame of (time) reference.

Here is a before and after plot. The before plot is terrible! (But flies fine.)  The after plot isn’t perfect, but might be about as good as it gets on a linux system.  Notice the difference in Y-scale between the two plots.  If you think your system is better than mine, log your IMU data at 100hz and plot the dt between samples and see for yourself.  In the following plots, the Y axis is dt time in seconds.  The X axis is elapsed run time in seconds.

imu_dt_beforedt using host’s timestamp when imu packet received.

imu_dt_afterdt when using a hybrid of host and apm2 timestamps.

Even with this fix, I see the host system’s main loop timing vary between 0.008 and 0.012 seconds per frame, occasionally even worse (100hz should ideally equal exactly 0.010 seconds.)  This is now far better than the system was doing previously, and far, far better than the PX4 does … but still not perfect.  There is always more work to do!

Conclusions

These mistakes (when finally discoveded) all led to important improvements with the AuraUAS system: better accelerometer calibration, better gyro resolution, better time step consistency with no dropped frames.  Will it help airplanes get from point A to point B more smoothly and more precisely?  Probably not in any externally visible way.  Mistakes?  I still make them 1000’s of times a day.  Lurking hidden mistakes?  Yes, those too.  My hope is that no matter what stage of life I find myself in, I’m always working for improvements, always vigilant to spot issues, and always focused on addressing issues when they are discovered.

Automated Movie Frame Extracting and Geotagging

This is a short tutorial on an automated method to extract and geotag movie frames.  One specific use case is that you have just flown a survey with your quad copter using a 2-axis gimbal pointing straight down, and a gopro action cam in movie mode.  Now you’d like to create a stitched map from your data using tools like pix4d or agisoft.

The most interesting part of this article is the method I have developed to correlate the frame timing of a movie with the aircraft’s flight data log.  This correlation process yields a result such that for any and every frame of the movie, I can find the exact corresponding time in the flight log, and for any time in the flight log, I can find the corresponding video frame.  Once this relationship is established, it is a simple matter to walk though the flight log and pull frames based on the desired conditions (for example, grab a frame at some time interval, while above some altitude AGL, and only when oriented +/- 10 degrees from north or south.)

Video Analysis

The first step of the process is to analyze the frame to frame motion in the video.

features

Example of feature detection in a single frame.

  1. For each video frame we run a feature detection algorithm (such as SIFT, SURF, Orb, etc.) and then compute the descriptors for each found feature.
  2. Find all the feature matches between frame “n-1” and frame “n”.  This is done using standard FLANN matching, followed by a RANSAC based homography matrix solver, and then discarding outliers.  This approach has a natural advantage of being able to ignore extraneous features from the prop or the nose of the aircraft because those features don’t fit into the overall consensus of the homography matrix.
  3. Given the set of matching features between frame “n-1” and frame “n”, I then compute a best fit rigid affine matrix transformation from the features locations (u, v) from one frame to the next.  The affine transformation can be decomposed into a rotational component, a translation (x, y) component, and a scale component.
  4. Finally I log the frame number, frame time (starting at t=0.0 for the first frame), and the rotation (deg), x translation (pixels), and y translation (pixels.)

The cool, tricky observation

I haven’t seen anyone else do anything like this before, so I’ll pretend I’ve come up with something new and cool.  I know there is never anything new under the sun, but it’s fun to rediscover things for oneself.

Use case #1: Iris quad copter with a two axis Tarot gimbal, and a go-pro hero 3 pointing straight down.  Because the gimbal is only 2-axis, the camera tracks the yaw motion of the quad copter exactly.  The camera is looking straight down, so camera roll rate should exactly match the quad copter’s yaw rate.  I have shown I can compute the frame to frame roll rate using computer vision techniques, and we can save the iris flight log.  If these two signal channels aren’t too noisy or biased relative to each other, perhaps we can find a way to correlate them and figure out the time offset.

iris-tarot

3DR Iris + Tarot 2-axis gimbal

Use case #2:  A Senior Telemaster fixed wing aircraft with a mobius action cam fixed to the airframe looking straight forward.  In this example, camera roll should exactly correlate to aircraft roll.  Camera x translation should map to aircraft yaw, and camera y translation should map to aircraft pitch.

IMG_20150804_073522

Senior Telemaster with forward looking camera.

In all cases this method requires that at least one of the camera axes is fixed relative to at least one of the aircraft axes.  If you are running a 3 axis gimbal you are out of luck … but perhaps with this method in mind and a bit of ingenuity alternative methods could be devised to find matching points in the video versus the flight log.

Flight data correlation

This is the easy part.  After processing the movie file, we  now have a log of the frame to frame motion.  We also have the flight log from the aircraft.  Here are the steps to correlate the two data logs.

overlay

Correlated sensor data streams.

  1. Load both data logs (movie log and flight log) into a big array and resample the data at a consistent interval.  I have found that resampling at 30hz seems to work well enough.  I have experimented with fitting a spline curve through lower rate data to smooth it out.  It makes the plots look prettier, but I’m sure does not improve the accuracy of the correlation.
  2. I coded this process up in python.  Luckily python (numpy) has a function that takes two time sequences as input and does brute force correlation.  It slides the one data stream forward against the other data stream and computes a correlation value for every possibly overlap.  This is why it is important to resample both data streams at the same fixed sample rate.
    ycorr = np.correlate(movie_interp[:,1], flight_interp[:,1], mode='full')
  3. When you plot out “ycorr”, you will hopefully see a spike in the plot and that should correspond to the best fit of the two data streams.

correlation

Plot of data overlay position vs. correlation.

Geotagging move frames

 

GOPR0002-009054

Raw Go-pro frame grab showing significant lens distortion.  Latitude = 44.69231071, Longitude = -93.06131655, Altitude = 322.1578

The important result of the correlation step is that we have now determined the exact offset in seconds between the movie log and the flight log.  We can use this to easily map a point in one data file to a point in the other data file.

Movie encoding formats are sequential and the compression algorithms require previous frames to generate the next frame.  Thus the geotagging script steps through the movie file frame by frame and finds the point in the flight log data file that matches.

For each frame that matches the extraction conditions, it is a simple matter to lookup the corresponding longitude, latitude, and altitude from the flight log.  My script provides an example of selecting movie frames based on conditions in the flight log.  I know that the flight was planned so the transacts were flown North/South and the target altitude was about 40m AGL.  I specifically coded the script to extract movie frames at a specified interval in seconds, but only consider frames taken when the quad copter was above 35m AGL and oriented within +/- 10 degrees of either North or South.  The script is written in python so it could easily be adjusted for other constraints.

The script writes each selected frame to disk using the opencv imwrite() function, and then uses the python “pyexiv2” module to write the geotag information into the exif header for that image.

geolocated

A screen grab from Pix4d showing the physical location of all the captured Go-pro movie frames.

Applications?

Aerial surveying and mapping

The initial use case for this code was to automate the process of extracting frames from a go-pro movie and geotagging them in preparation for handing the image set over to pix4d for stitching and mapping.

gopro-stitch

Final stitch result from 120 geotagged gopro movie frames.

Using video as a truth reference to analyze sensor quality

It is interesting to see how accurately the video roll rate corresponds to the IMU gyro roll rate (assume a forward looking camera now.)  It is also interesting in plots to see how the two data streams track exactly for some periods of time, but diverge by some slowly varying bias for other periods of time.  I believe this shows the variable bias of MEMS gyro sensors.  It would be interesting to run down this path a little further and see if the bias correlates to g force in a coupled axis?

Visual odometry and real time mapping

Given feature detection and matching from one frame to the next, knowledge of the camera pose at each frame, opencv pnp() and triangulation() functions, and a working bundle adjuster … what could be done to map the surface or compute visual odometry during a gps outage?

Source Code

The source code for all my image analysis experimentation can be found at the University of Minnesota UAV Lab github page.  It is distributed under the MIT open-source license:

https://github.com/UASLab/ImageAnalysis

Comments or questions?

I’d love to see your comments or questions in the comments section at the end of this page!

Image Stitching Tutorial Part #2: Direct Georeferencing

What is Direct Georeferencing?

uav-for-PA_01-image-collection

The process of mapping a pixel in an image to a real world latitude, longitude, and altitude is called georeferencing.  When UAS flight data is tightly integrated with the camera and imagery data, the aerial imagery can be directly placed and oriented on a map.   Any image feature can be directly located.  All of this can be done without needing to feature detect, feature match, and image stitch.  The promise of “direct georeferencing” is the ability to provide useful maps and actionable data immediately after the conclusion of a flight.

Typically a collection of overlapping images (possibly tagged with the location of the camera when the image was snapped) are stitched together to form a mosaic.  Ground control points can be surveyed and located in the mosaic to georectify the entire image set.   However the process of stitching imagery together can be very time consuming and specific images sets can fail to stitch well for a variety of reasons.

If a good estimate for the location and orientation of the camera is available for each picture, and if a good estimate of the ground surface is available, the tedious process of image stitching can be skipped and image features can be georeferenced directly.

Where is Direct Georeferencing Suitable for use?

stressed_trees

There are a number of reasons to skip the traditional image stitching pipeline and use direct georeferencing.  I list a few here, but if the readers can think of other use cases, I would love to hear your thoughts and expand this list!

  • Any situation where we are willing to trade a seamless ‘nice’ looking mosaic for fast results.  i.e. a farmer in a field that would like to make a same-day decision rather than wait for a day or two for image stitching results to be computed.
  • Surveys or counting.  One use case I have worked on is marine surveys.  In this use case the imagery is collected over open-ocean with no stable features to stitch.  Instead we were more interested on finding things that were not water and getting a quick but accurate location estimate for them.  Farmers might want to count livestock, land managers might want to locate dead trees, researchers might want to count bird populations on a remote island.

debris1

debris

How can Direct Georeferencing Improve the Traditional Image Stitching Pipeline?

There are some existing commercial image stitching applications that are very good at what they do.  However, they are closed-source and don’t give up their secrets easily.  Thus it is hard to do an apples-to-apples comparison with commercial tools to evaluate how (and how much) direct georeferencing can improve the traditional image stitching pipeline.  With that in mind I will forge ahead and suggest several ways I believe direct georeferencing can improve the traditional methods:

  • Direct georeferencing provides an initial 3d world coordinate estimate for every detected feature before any matching or stitching work is performed.
  • The 3d world coordinates of features can be used to pre-filter match sets between images.  When doing an n vs. n image compare to find matching image pairs, we can compare only images with feature sets that overlap in world coordinates, and then only compare the overlapping subset of features.  This speeds up the feature matching process by reducing the number of individual feature comparisons.  This increases the robustness of the feature matching process by reducing the number of potential similar features in the comparison set.  And this helps find potential match pairs that other applications might miss.
  • After an initial set of potential feature matches is found between a pair of images, these features must be further evaluated and filtered to remove false matches.  There are a number of approaches for filtering out incorrect matches, but with direct georeferencing we can add real world proximity to our set of strategies for eliminating bad matches.
  • Once the entire match set is computed for the entire image set, the 3d world coordinate for each matched feature can be further refined by averaging the estimate from each matching image together.
  • When submitting point and camera data to a bundle adjustment algorithm, we can provide our positions already in 3d world coordinates.  We don’t need to build up an initial estimate in some arbitrary camera coordinate system where each image’s features are positioned relative to neighbors.  Instead we can start with a good initial guess for the location of all our features.
  • When the bundle adjustment algorithm finishes.  We can compare the new point location estimates against the original estimates and look for features that have moved implausibly far.  This could be evidence of remaining outliers.

Throughout the image stitching pipeline, it is critical to create a set of image feature matches that link all the images together and cover the overlapping portions of each image pair as fully as possible.  False matches can cause ugly imperfections in the final stitched result and they can cause the bundle adjustment algorithm to fail so it is critical to find a good set of feature matches.  Direct georeferencing can improve the critical feature matching process.

What Information is Required for Direct Georeferencing?

  • An accurate camera position and orientation for each image.  This may require the flight data log from your autopilot and an accurate (sub second) time stamp when the image was snapped.
  • An estimate of the ground elevation or terrain model.
  • Knowledge of your camera’s lens calibration (“K” matrix.)  This encompasses the field of view of your lens (sensor dimensions and focal length) as well as the sensor width and height in pixels.
  • Knowledge of your camera’s lens distortion parameters.  Action cameras like a gopro or mobius have significant lens distortion that must be accounted for.

Equipment to Avoid

I don’t intend this section to be negative, every tool has strengths and times when it is a good choice.  However, it is important to understand what equipment choices works better and what choices may present challenges.

  • Any camera that makes it difficult to snap a picture precisely (< 0.1 seconds) from the trigger request.  Autofocus commercial cameras can introduce random latency between the camera trigger and the actual photo.
  •  Rolling shutter cameras.  This is just about every commercial off the shelf camera sadly, but rolling shutter introduces warping into the image which can add uncertainty to the results.  This can be partially mitigated by setting your camera to a very high shutter speed (i.e. 1/800 or 1/1000th of a second.)
  • Cameras with slow shutter speeds or cameras that do not allow you to set your shutter speed or other parameters.
  • Any camera mounted to an independent gimbal.  A gimbal works nice for stable video, but if it separates the camera orientation from the aircraft orientation, then we can no longer use aircraft orientation to compute camera orientation.
  • Any flight computer that doesn’t let you download a complete flight log that includes real world time stamp, latitude, longitude, altitude, roll, pitch, and yaw.

The important point I am attempting to communicate is that tight integration between the camera and the flight computer is an important aspect for direct georeferencing.  Strapping a gimbaled action cam to a commercial quad-copter very likely will not allow you to extract all the information required for direct georeferencing.

 

Image Stitching Tutorial Part #1: Introduction

features

Background

During the summer of 2014 I began investigating image stitching techniques and technologies for a NOAA sponsored UAS marine survey project.  In the summer of 2015 I was hired by the University of Minnesota Department of Aerospace Engineering and Mechanics to work on a Precision Agriculture project that also involves UAS’s and aerial image stitching.

Over the past few months I have developed a functional open-source image stitching pipeline written in python and opencv.  It is my intention with this series of blog postings to introduce this work and further explain our approach to aerial image processing and stitching.

Any software development project is a journey of discovery and education, so I would love to hear your thoughts, feedback, and questions in the comments area of any of these posts.  The python code described here will be released under the MIT open-source software license (one of my to-do list items is to publish this project code, so that will happen “soon.”)
feature-matching

 Why?

The world already has several high quality commercial image stitching tools as well as several cloud based systems that are free to use.  Why develop yet another image stitching pipeline?  There are several reasons we began putting this software tool chain together.

  • We are looking at the usefulness and benefits of ‘direct georeferencing.’  If we have accurate camera pose information (time, location, and camera orientation of each image), then how can this improve the image analysis, feature detection, feature matching, stitching, and rendering process?
  • One of the the core strengths of the UMN Aerospace Engineering Department is a high quality 15-state kalman filter attitude determination system.  This system uses inertial sensors (gyros and accelerometers) in combination with a GPS to accurately estimate an aircraft’s ‘true’ orientation and position.  Our department is uniquely positioned to provide a high quality camera pose estimate and thus examine ‘direct georeferencing’ based image processing.
  • Commercial software and closed source cloud solutions do not enable the research community to easily ask questions and test ideas and theories.
  • We hope to quantify the sensor quality required to perform useful direct georeferencing as well as the various sources of uncertainty that can influence the results.
  • We would love to involve the larger community in our work, and we anticipate there will be some wider interest in free/open image processing and stitching tools that anyone can modify and run on their own computer.

Outline

I will be posting new tutorials in this series as they are written.  Here is a quick look ahead at what topics I plan to cover:

  • Direct Georeferencing
  • Image stitching basics
  • Introduction to the open-source software tool chain
  • Aircraft vs. camera poses and directly visualizing your georeferenced data set
  • Feature detection
  • Feature matching
  • Sparse Bundle Adjustment
  • Seamless 3D and 2D image mosaics, DEM’s, Triangle meshes, etc.

Throughout the image collection and image stitching process there is art, science, engineering, math, software, hardware, aircraft, skill, and a maybe bit of luck once in a while (!) that all come together in order to produce a successful aerial imaging result.

Software Download

The software referenced in this tutorial series is licensed with the MIT license and available on the University of Minnestoa UAV Lab public github page under the ImageAnalysis repository.

Credits

Flight Milestones

Congratulations!

Congrats ATI Resolution 3, Hobby Lobby Senior Telemaster, Hobbyking Skywalker, and Avior-lite autopilot on your recent milestones!

IMG_20130923_183033 IMG_20150804_073522 IMG_20150804_093541 IMG_20150801_135249

Avior-lite (beaglebone + apm2 hybrid) autopilot:

  • 300th logged flight
  • 7000+ logged flight minutes (117.8 hours)
  • 6400+ fully autonomous flight minutes (107.2 hours)
  • 2895 nautical miles flown (3332 miles, 5362 km)

Hobby Lobby Senior Telemaster (8′ wing span)

  • Actively flight testing autopilot hardware and software changes since 2007!
  • 200th logged flight.
  • 5013 logged flight minutes (83.5 hours)
  • 4724 fully autonomous flight minutes (78.7 hours)
  • 2015 nautical miles flown (2319 miles, 3733 km)

Today (October 7, 2015) I logged the 300th avior-lite flight and simultaneously logged the 200th flight on my venerable Senior Telemaster.  I realize these are just numbers, and they wouldn’t be big numbers for a full scale aircraft operation or even a production uav operation, but it represents my personal effort in the UAV field.

I’m proud of a mishap-free 2015 flying season so far!  (Ok, err, well one mishap on the very first launch of the skywalker … grrr … pilot error … and fixable thankfully.)

Enjoy the fall colors and keep flying!

fall-colors-skywalker