Thursday, January 13, 2011

Fast Unwrap

On towards my first Matlab entry:

I often have to 'unwrap' longitudinal data.
I have a vector a matrix of longitudinal data
lon(m,n)
now these get displaced, shifted or what so ever and the original matrix has elements that are >180 or <-180. (Or if you prefer the 0 to 360 space: <0 or >360, or pi respectively)

To plot these data I usually have to unwrap them which I did with:


lon(lon>180) = lon(lon>180) -360;
lon(lon<=-180) = lon(lon<=-180) +360;


but doing so several times, hundreds of times takes some time. Following approach does not need to compare sizes and is significant faster:


lon = lon - 360.*sign(lon).*(1+sign(abs(lon+eps)-180))./2;


the +eps is necessary to prevent sign(0), assuming data is not close to machine precision.
There is no need to add an eps in the middle of the 2nd term, since the last part of term two is zero in case lon is close to zero.
Now you can improve this code to the following, which actually is substantially faster, for the fact that two scalars get added prior to the addition of a scalar to a vector, compared to the above case in which twice a scalar is added to a vector, but with the limitation that you now only check for < and > at then ends or  => and <= at both ends, no longer only the = in one direction:


lon = lon - 360.*sign(lon).*(1+sign(abs(lon)-(180-eps)))./2;

In case you are looking at (lon>=180) and (lon<=-180) data, just turn the -eps into a +eps.

But even in case you know your longitudinal data is displaced only into one direction and your code looks like:


lon(lon>180) = lon(lon>180) -360;


the code above does 'extends' to:


lon = lon - 360.*(1+sign(lon-(180-eps)))./2;


and similar approach for lon >= 180 instead just > as above: change the sign in front of the eps.
Still this code can be excuted much faster if you rewrite it to:


lon = lon - 180 + 180*sign(lon-(180-eps)));

Just how much faster does this code execute, well on my machine you need about 33% less computation time for data with larger and smaller values with the new approach and
for the case with data only displaced into one direction it is about 25%.
But with data on a sphere, cylinder or circle where you do unwrap regularly it can accumulate to a significant amount.

I didn't come up with an even faster way - maybe someone else has an even more sophisticated idea.

oh - I forgot, in case you wondered how this compares to the Matlab internal 'unwrap' - it is way way WAY faster!

No comments:

Post a Comment