Saturday, March 19, 2016

Error in NASA code

NASA World Wind java code is available on the internet. While looking at the their SDK2 EGM96.java file, I discovered that the simple two dimensional linear (bilenear) interpolation code has an error. The original code:
...
double ul = this.gePostOffset(topRow, leftCol);
double ll = this.gePostOffset(bottomRow, leftCol);
double lr = this.gePostOffset(bottomRow, rightCol);
double ur = this.gePostOffset(topRow, rightCol);

double u = (lon - lonLeft) / INTERVAL.degrees;
double v = (latTop - lat) / INTERVAL.degrees;

double pll = (1.0 - u) * (1.0 - v);
double plr = u * (1.0 - v);
double pur = u * v;
double pul = (1.0 - u) * v;

double offset = pll * ll + plr * lr + pur * ur + pul * ul;
...

I prepared a schematic to visualize the algorithm:


As you can see at the end of above code snippet, the offset formula starts with pll*ll. pll is (1.0 - u) * (1.0 - v). The correct multiplier of the ll term has to be v*(1-u).

A quick sanity check: Consider the case when lat=latBottom and lon=lonLeft, i.e. the lower left corner. In that case we expect the offset to be equal to ll. Plugging in the values we get v=1, 1-v = 0 and u=0, 1-u = 1. If we use the original (wrong) formula we get offset = 1*0*ll + 0*0*lr + 0*1*ur + 1*0*lu = 0 which is clearly wrong. If we use the corrected formula we get offset = 1*1*ll + 0*0*lr + 0*1*ur + 1*0*lu = ll which is the expected result.

The easiest way to fix the code is to swap the topRow and bottomRow indices when calling getPostOffset() function.

I guess the error was not detected since results do not differ too much because, although the points are used in wrong order, they are close to each other nonetheless, so the result might not have looked suspicious.

I tried to log a bug report on their Jira site, but I could not since I don't have an account. I was able to file a bug report on GitHub.

Lessons learnt:
  • Thoroughly test 3rd party code, even if it is from NASA or Mathworks. Every code is guilty unless proven innocent. Ignore this advice and you will find yourself chasing strange errors for a very long time.
  • Separate interpolation code into a function, write lots of unit test for that function to verify that it works correctly.

3 comments:

Anonymous said...

Can you explain the offset formula? It does't seem bilinear interpolation formula. With bilinear interpolation I get different results.

Samil Korkmaz said...

Looking at the Wikipedia article for bilinear interpolation, we can map our problem here as follows: ll = Q11, lr = Q21, ur = Q22, ul = Q12. 1-u = (x2-x)/(x2-x1), v = (y2-y)/(y2-y1) and so on. The first term v*(1-u)*ll will be (y2-y)/(y2-y1)*(x2-x)/(x2-x1)*Q11 which is the same first term in Wikipedia. If you replace the rest, you will see that the formula here and the formula in Wikipedia are the same. Please check your calculations.

Zach Glueckert said...

Samil, thanks for spotting this issue with NASA WorldWindJava! Your analysis was correct and we've since updated WorldWindJava to incorporate a fix. Check out https://github.com/NASAWorldWind/WorldWindJava and https://github.com/NASAWorldWind/ for other versions of WorldWind and the latest code updates. Thanks again!