src/viktrack.c

   1 /*
   2  * viking -- GPS Data and Topo Analyzer, Explorer, and Manager
   3  *
   4  * Copyright (C) 2003-2005, Evan Battaglia <gtoevan@gmx.net>
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 2 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19  *
  20  */
  21 #ifdef HAVE_CONFIG_H
  22 #include "config.h"
  23 #endif
  24
  25 #include <glib.h>
  26 #include <time.h>
  27 #include <stdlib.h>
  28 #ifdef HAVE_STRING_H
  29 #include <string.h>
  30 #endif
  31 #ifdef HAVE_MATH_H
  32 #include <math.h>
  33 #endif
  34
  35 #include "coords.h"
  36 #include "vikcoord.h"
  37 #include "viktrack.h"
  38 #include "globals.h"
  39 #include "dems.h"
  40
  41 VikTrack *vik_track_new()
  42 {
  43   VikTrack *tr = g_malloc0 ( sizeof ( VikTrack ) );
  44   tr->ref_count = 1;
  45   return tr;
  46 }
  47
  48 void vik_track_set_comment_no_copy(VikTrack *tr, gchar *comment)
  49 {
  50   if ( tr->comment )
  51     g_free ( tr->comment );
  52   tr->comment = comment;
  53 }
  54
  55
  56 void vik_track_set_name(VikTrack *tr, const gchar *name)
  57 {
  58   if ( tr->name )
  59     g_free ( tr->name );
  60
  61   if ( name && name[0] != '\0' )
  62     tr->name = g_strdup(name);
  63   else
  64     tr->name = NULL;
  65 }
  66
  67 void vik_track_set_comment(VikTrack *tr, const gchar *comment)
  68 {
  69   if ( tr->comment )
  70     g_free ( tr->comment );
  71
  72   if ( comment && comment[0] != '\0' )
  73     tr->comment = g_strdup(comment);
  74   else
  75     tr->comment = NULL;
  76 }
  77
  78 void vik_track_ref(VikTrack *tr)
  79 {
  80   tr->ref_count++;
  81 }
  82
  83 void vik_track_set_property_dialog(VikTrack *tr, GtkWidget *dialog)
  84 {
  85   /* Warning: does not check for existing dialog */
  86   tr->property_dialog = dialog;
  87 }
  88
  89 void vik_track_clear_property_dialog(VikTrack *tr)
  90 {
  91   tr->property_dialog = NULL;
  92 }
  93
  94 void vik_track_free(VikTrack *tr)
  95 {
  96   if ( tr->ref_count-- > 1 )
  97     return;
  98
  99   if ( tr->name )
 100     g_free ( tr->name );
 101   if ( tr->comment )
 102     g_free ( tr->comment );
 103   g_list_foreach ( tr->trackpoints, (GFunc) g_free, NULL );
 104   g_list_free( tr->trackpoints );
 105   if (tr->property_dialog)
 106     if ( GTK_IS_WIDGET(tr->property_dialog) )
 107       gtk_widget_destroy ( GTK_WIDGET(tr->property_dialog) );
 108   g_free ( tr );
 109 }
 110
 111 VikTrack *vik_track_copy ( const VikTrack *tr )
 112 {
 113   VikTrack *new_tr = vik_track_new();
 114   VikTrackpoint *new_tp;
 115   GList *tp_iter = tr->trackpoints;
 116   new_tr->visible = tr->visible;
 117   new_tr->trackpoints = NULL;
 118   while ( tp_iter )
 119   {
 120     new_tp = g_malloc ( sizeof ( VikTrackpoint ) );
 121     *new_tp = *((VikTrackpoint *)(tp_iter->data));
 122     new_tr->trackpoints = g_list_append ( new_tr->trackpoints, new_tp );
 123     tp_iter = tp_iter->next;
 124   }
 125   vik_track_set_name(new_tr,tr->name);
 126   vik_track_set_comment(new_tr,tr->comment);
 127   return new_tr;
 128 }
 129
 130 VikTrackpoint *vik_trackpoint_new()
 131 {
 132   VikTrackpoint *tp = g_malloc0(sizeof(VikTrackpoint));
 133   tp->speed = NAN;
 134   tp->course = NAN;
 135   tp->altitude = VIK_DEFAULT_ALTITUDE;
 136   tp->hdop = VIK_DEFAULT_DOP;
 137   tp->vdop = VIK_DEFAULT_DOP;
 138   tp->pdop = VIK_DEFAULT_DOP;
 139   return tp;
 140 }
 141
 142 void vik_trackpoint_free(VikTrackpoint *tp)
 143 {
 144   g_free(tp);
 145 }
 146
 147 VikTrackpoint *vik_trackpoint_copy(VikTrackpoint *tp)
 148 {
 149   VikTrackpoint *rv = vik_trackpoint_new();
 150   *rv = *tp;
 151   return rv;
 152 }
 153
 154 gdouble vik_track_get_length(const VikTrack *tr)
 155 {
 156   gdouble len = 0.0;
 157   if ( tr->trackpoints )
 158   {
 159     GList *iter = tr->trackpoints->next;
 160     while (iter)
 161     {
 162       if ( ! VIK_TRACKPOINT(iter->data)->newsegment )
 163         len += vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 164                                 &(VIK_TRACKPOINT(iter->prev->data)->coord) );
 165       iter = iter->next;
 166     }
 167   }
 168   return len;
 169 }
 170
 171 gdouble vik_track_get_length_including_gaps(const VikTrack *tr)
 172 {
 173   gdouble len = 0.0;
 174   if ( tr->trackpoints )
 175   {
 176     GList *iter = tr->trackpoints->next;
 177     while (iter)
 178     {
 179       len += vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 180                               &(VIK_TRACKPOINT(iter->prev->data)->coord) );
 181       iter = iter->next;
 182     }
 183   }
 184   return len;
 185 }
 186
 187 gulong vik_track_get_tp_count(const VikTrack *tr)
 188 {
 189   gulong num = 0;
 190   GList *iter = tr->trackpoints;
 191   while ( iter )
 192   {
 193     num++;
 194     iter = iter->next;
 195   }
 196   return num;
 197 }
 198
 199 gulong vik_track_get_dup_point_count ( const VikTrack *tr )
 200 {
 201   gulong num = 0;
 202   GList *iter = tr->trackpoints;
 203   while ( iter )
 204   {
 205     if ( iter->next && vik_coord_equals ( &(VIK_TRACKPOINT(iter->data)->coord),
 206                        &(VIK_TRACKPOINT(iter->next->data)->coord) ) )
 207       num++;
 208     iter = iter->next;
 209   }
 210   return num;
 211 }
 212
 213 void vik_track_remove_dup_points ( VikTrack *tr )
 214 {
 215   GList *iter = tr->trackpoints;
 216   while ( iter )
 217   {
 218     if ( iter->next && vik_coord_equals ( &(VIK_TRACKPOINT(iter->data)->coord),
 219                        &(VIK_TRACKPOINT(iter->next->data)->coord) ) )
 220     {
 221       g_free ( iter->next->data );
 222       tr->trackpoints = g_list_delete_link ( tr->trackpoints, iter->next );
 223     }
 224     else
 225       iter = iter->next;
 226   }
 227 }
 228
 229 guint vik_track_get_segment_count(const VikTrack *tr)
 230 {
 231   guint num = 1;
 232   GList *iter = tr->trackpoints;
 233   if ( !iter )
 234     return 0;
 235   while ( (iter = iter->next) )
 236   {
 237     if ( VIK_TRACKPOINT(iter->data)->newsegment )
 238       num++;
 239   }
 240   return num;
 241 }
 242
 243 VikTrack **vik_track_split_into_segments(VikTrack *t, guint *ret_len)
 244 {
 245   VikTrack **rv;
 246   VikTrack *tr;
 247   guint i;
 248   guint segs = vik_track_get_segment_count(t);
 249   GList *iter;
 250
 251   if ( segs < 2 )
 252   {
 253     *ret_len = 0;
 254     return NULL;
 255   }
 256
 257   rv = g_malloc ( segs * sizeof(VikTrack *) );
 258   tr = vik_track_copy ( t );
 259   rv[0] = tr;
 260   iter = tr->trackpoints;
 261
 262   i = 1;
 263   while ( (iter = iter->next) )
 264   {
 265     if ( VIK_TRACKPOINT(iter->data)->newsegment )
 266     {
 267       iter->prev->next = NULL;
 268       iter->prev = NULL;
 269       rv[i] = vik_track_new();
 270       // TODO: consider new naming strategy here
 271       if ( tr->name )
 272         vik_track_set_name ( rv[i], tr->name );
 273       if ( tr->comment )
 274         vik_track_set_comment ( rv[i], tr->comment );
 275       rv[i]->visible = tr->visible;
 276       rv[i]->trackpoints = iter;
 277       i++;
 278     }
 279   }
 280   *ret_len = segs;
 281   return rv;
 282 }
 283
 284 void vik_track_reverse ( VikTrack *tr )
 285 {
 286   GList *iter;
 287   tr->trackpoints = g_list_reverse(tr->trackpoints);
 288
 289   /* fix 'newsegment' */
 290   iter = g_list_last ( tr->trackpoints );
 291   while ( iter )
 292   {
 293     if ( ! iter->next ) /* last segment, was first, cancel newsegment. */
 294       VIK_TRACKPOINT(iter->data)->newsegment = FALSE;
 295     if ( ! iter->prev ) /* first segment by convention has newsegment flag. */
 296       VIK_TRACKPOINT(iter->data)->newsegment = TRUE;
 297     else if ( VIK_TRACKPOINT(iter->data)->newsegment && iter->next )
 298     {
 299       VIK_TRACKPOINT(iter->next->data)->newsegment = TRUE;
 300       VIK_TRACKPOINT(iter->data)->newsegment = FALSE;
 301     }
 302     iter = iter->prev;
 303   }
 304 }
 305
 306 gdouble vik_track_get_average_speed(const VikTrack *tr)
 307 {
 308   gdouble len = 0.0;
 309   guint32 time = 0;
 310   if ( tr->trackpoints )
 311   {
 312     GList *iter = tr->trackpoints->next;
 313     while (iter)
 314     {
 315       if ( VIK_TRACKPOINT(iter->data)->has_timestamp &&
 316           VIK_TRACKPOINT(iter->prev->data)->has_timestamp &&
 317           (! VIK_TRACKPOINT(iter->data)->newsegment) )
 318       {
 319         len += vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 320                                 &(VIK_TRACKPOINT(iter->prev->data)->coord) );
 321         time += ABS(VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(iter->prev->data)->timestamp);
 322       }
 323       iter = iter->next;
 324     }
 325   }
 326   return (time == 0) ? 0 : ABS(len/time);
 327 }
 328
 329 /**
 330  * Based on a simple average speed, but with a twist - to give a moving average.
 331  *  . GPSs often report a moving average in their statistics output
 332  *  . bicycle speedos often don't factor in time when stopped - hence reporting a moving average for speed
 333  *
 334  * Often GPS track will record every second but not when stationary
 335  * This method doesn't use samples that differ over the specified time limit - effectively skipping that time chunk from the total time
 336  *
 337  * Suggest to use 60 seconds as the stop length (as the default used in the TrackWaypoint draw stops factor)
 338  */
 339 gdouble vik_track_get_average_speed_moving (const VikTrack *tr, int stop_length_seconds)
 340 {
 341   gdouble len = 0.0;
 342   guint32 time = 0;
 343   if ( tr->trackpoints )
 344   {
 345     GList *iter = tr->trackpoints->next;
 346     while (iter)
 347     {
 348       if ( VIK_TRACKPOINT(iter->data)->has_timestamp &&
 349           VIK_TRACKPOINT(iter->prev->data)->has_timestamp &&
 350           (! VIK_TRACKPOINT(iter->data)->newsegment) )
 351       {
 352         if ( ( VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(iter->prev->data)->timestamp ) < stop_length_seconds ) {
 353           len += vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 354                                   &(VIK_TRACKPOINT(iter->prev->data)->coord) );
 355
 356           time += ABS(VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(iter->prev->data)->timestamp);
 357         }
 358       }
 359       iter = iter->next;
 360     }
 361   }
 362   return (time == 0) ? 0 : ABS(len/time);
 363 }
 364
 365 gdouble vik_track_get_max_speed(const VikTrack *tr)
 366 {
 367   gdouble maxspeed = 0.0, speed = 0.0;
 368   if ( tr->trackpoints )
 369   {
 370     GList *iter = tr->trackpoints->next;
 371     while (iter)
 372     {
 373       if ( VIK_TRACKPOINT(iter->data)->has_timestamp &&
 374           VIK_TRACKPOINT(iter->prev->data)->has_timestamp &&
 375           (! VIK_TRACKPOINT(iter->data)->newsegment) )
 376       {
 377         speed =  vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord), &(VIK_TRACKPOINT(iter->prev->data)->coord) )
 378                  / ABS(VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(iter->prev->data)->timestamp);
 379         if ( speed > maxspeed )
 380           maxspeed = speed;
 381       }
 382       iter = iter->next;
 383     }
 384   }
 385   return maxspeed;
 386 }
 387
 388 void vik_track_convert ( VikTrack *tr, VikCoordMode dest_mode )
 389 {
 390   GList *iter = tr->trackpoints;
 391   while (iter)
 392   {
 393     vik_coord_convert ( &(VIK_TRACKPOINT(iter->data)->coord), dest_mode );
 394     iter = iter->next;
 395   }
 396 }
 397
 398 /* I understood this when I wrote it ... maybe ... Basically it eats up the
 399  * proper amounts of length on the track and averages elevation over that. */
 400 gdouble *vik_track_make_elevation_map ( const VikTrack *tr, guint16 num_chunks )
 401 {
 402   gdouble *pts;
 403   gdouble total_length, chunk_length, current_dist, current_area_under_curve, current_seg_length, dist_along_seg = 0.0;
 404   gdouble altitude1, altitude2;
 405   guint16 current_chunk;
 406   gboolean ignore_it = FALSE;
 407
 408   GList *iter = tr->trackpoints;
 409
 410   if (!iter || !iter->next) /* zero- or one-point track */
 411           return NULL;
 412
 413   { /* test if there's anything worth calculating */
 414     gboolean okay = FALSE;
 415     while ( iter )
 416     {
 417       if ( VIK_TRACKPOINT(iter->data)->altitude != VIK_DEFAULT_ALTITUDE ) {
 418         okay = TRUE; break;
 419       }
 420       iter = iter->next;
 421     }
 422     if ( ! okay )
 423       return NULL;
 424   }
 425
 426   iter = tr->trackpoints;
 427
 428   g_assert ( num_chunks < 16000 );
 429
 430   pts = g_malloc ( sizeof(gdouble) * num_chunks );
 431
 432   total_length = vik_track_get_length_including_gaps ( tr );
 433   chunk_length = total_length / num_chunks;
 434
 435   /* Zero chunk_length (eg, track of 2 tp with the same loc) will cause crash */
 436   if (chunk_length <= 0) {
 437     g_free(pts);
 438     return NULL;
 439   }
 440
 441   current_dist = 0.0;
 442   current_area_under_curve = 0;
 443   current_chunk = 0;
 444   current_seg_length = 0;
 445
 446   current_seg_length = vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 447       &(VIK_TRACKPOINT(iter->next->data)->coord) );
 448   altitude1 = VIK_TRACKPOINT(iter->data)->altitude;
 449   altitude2 = VIK_TRACKPOINT(iter->next->data)->altitude;
 450   dist_along_seg = 0;
 451
 452   while ( current_chunk < num_chunks ) {
 453
 454     /* go along current seg */
 455     if ( current_seg_length && (current_seg_length - dist_along_seg) > chunk_length ) {
 456       dist_along_seg += chunk_length;
 457
 458       /*        /
 459        *   pt2 *
 460        *      /x       altitude = alt_at_pt_1 + alt_at_pt_2 / 2 = altitude1 + slope * dist_value_of_pt_inbetween_pt1_and_pt2
 461        *     /xx   avg altitude = area under curve / chunk len
 462        *pt1 *xxx   avg altitude = altitude1 + (altitude2-altitude1)/(current_seg_length)*(dist_along_seg + (chunk_len/2))
 463        *   / xxx
 464        *  /  xxx
 465        **/
 466
 467       if ( ignore_it )
 468         // Seemly can't determine average for this section - so use last known good value (much better than just sticking in zero)
 469         pts[current_chunk] = altitude1;
 470       else
 471         pts[current_chunk] = altitude1 + (altitude2-altitude1)*((dist_along_seg - (chunk_length/2))/current_seg_length);
 472
 473       current_chunk++;
 474     } else {
 475       /* finish current seg */
 476       if ( current_seg_length ) {
 477         gdouble altitude_at_dist_along_seg = altitude1 + (altitude2-altitude1)/(current_seg_length)*dist_along_seg;
 478         current_dist = current_seg_length - dist_along_seg;
 479         current_area_under_curve = current_dist*(altitude_at_dist_along_seg + altitude2)*0.5;
 480       } else { current_dist = current_area_under_curve = 0; } /* should only happen if first current_seg_length == 0 */
 481
 482       /* get intervening segs */
 483       iter = iter->next;
 484       while ( iter && iter->next ) {
 485         current_seg_length = vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 486             &(VIK_TRACKPOINT(iter->next->data)->coord) );
 487         altitude1 = VIK_TRACKPOINT(iter->data)->altitude;
 488         altitude2 = VIK_TRACKPOINT(iter->next->data)->altitude;
 489         ignore_it = VIK_TRACKPOINT(iter->next->data)->newsegment;
 490
 491         if ( chunk_length - current_dist >= current_seg_length ) {
 492           current_dist += current_seg_length;
 493           current_area_under_curve += current_seg_length * (altitude1+altitude2) * 0.5;
 494           iter = iter->next;
 495         } else {
 496           break;
 497         }
 498       }
 499
 500       /* final seg */
 501       dist_along_seg = chunk_length - current_dist;
 502       if ( ignore_it || ( iter && !iter->next ) ) {
 503         pts[current_chunk] = current_area_under_curve / current_dist;
 504         if (!iter->next) {
 505           int i;
 506           for (i = current_chunk + 1; i < num_chunks; i++)
 507             pts[i] = pts[current_chunk];
 508           break;
 509         }
 510       }
 511       else {
 512         current_area_under_curve += dist_along_seg * (altitude1 + (altitude2 - altitude1)*dist_along_seg/current_seg_length);
 513         pts[current_chunk] = current_area_under_curve / chunk_length;
 514       }
 515
 516       current_dist = 0;
 517       current_chunk++;
 518     }
 519   }
 520
 521   return pts;
 522 }
 523
 524
 525 void vik_track_get_total_elevation_gain(const VikTrack *tr, gdouble *up, gdouble *down)
 526 {
 527   gdouble diff;
 528   *up = *down = 0;
 529   if ( tr->trackpoints && VIK_TRACKPOINT(tr->trackpoints->data)->altitude != VIK_DEFAULT_ALTITUDE )
 530   {
 531     GList *iter = tr->trackpoints->next;
 532     while (iter)
 533     {
 534       diff = VIK_TRACKPOINT(iter->data)->altitude - VIK_TRACKPOINT(iter->prev->data)->altitude;
 535       if ( diff > 0 )
 536         *up += diff;
 537       else
 538         *down -= diff;
 539       iter = iter->next;
 540     }
 541   } else
 542     *up = *down = VIK_DEFAULT_ALTITUDE;
 543 }
 544
 545 gdouble *vik_track_make_gradient_map ( const VikTrack *tr, guint16 num_chunks )
 546 {
 547   gdouble *pts;
 548   gdouble *altitudes;
 549   gdouble total_length, chunk_length, current_gradient;
 550   gdouble altitude1, altitude2;
 551   guint16 current_chunk;
 552
 553   g_assert ( num_chunks < 16000 );
 554
 555   total_length = vik_track_get_length_including_gaps ( tr );
 556   chunk_length = total_length / num_chunks;
 557
 558   /* Zero chunk_length (eg, track of 2 tp with the same loc) will cause crash */
 559   if (chunk_length <= 0) {
 560     return NULL;
 561   }
 562
 563   altitudes = vik_track_make_elevation_map (tr, num_chunks);
 564   if (altitudes == NULL) {
 565     return NULL;
 566   }
 567
 568   current_gradient = 0.0;
 569   pts = g_malloc ( sizeof(gdouble) * num_chunks );
 570   for (current_chunk = 0; current_chunk < (num_chunks - 1); current_chunk++) {
 571     altitude1 = altitudes[current_chunk];
 572     altitude2 = altitudes[current_chunk + 1];
 573     current_gradient = 100.0 * (altitude2 - altitude1) / chunk_length;
 574
 575     pts[current_chunk] = current_gradient;
 576   }
 577
 578   pts[current_chunk] = current_gradient;
 579
 580   return pts;
 581 }
 582
 583 /* by Alex Foobarian */
 584 gdouble *vik_track_make_speed_map ( const VikTrack *tr, guint16 num_chunks )
 585 {
 586   gdouble *v, *s, *t;
 587   gdouble duration, chunk_dur;
 588   time_t t1, t2;
 589   int i, pt_count, numpts, index;
 590   GList *iter;
 591
 592   if ( ! tr->trackpoints )
 593     return NULL;
 594
 595   g_assert ( num_chunks < 16000 );
 596
 597   t1 = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 598   t2 = VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->timestamp;
 599   duration = t2 - t1;
 600
 601   if ( !t1 || !t2 || !duration )
 602     return NULL;
 603
 604   if (duration < 0) {
 605     g_warning("negative duration: unsorted trackpoint timestamps?");
 606     return NULL;
 607   }
 608   pt_count = vik_track_get_tp_count(tr);
 609
 610   v = g_malloc ( sizeof(gdouble) * num_chunks );
 611   chunk_dur = duration / num_chunks;
 612
 613   s = g_malloc(sizeof(double) * pt_count);
 614   t = g_malloc(sizeof(double) * pt_count);
 615
 616   iter = tr->trackpoints->next;
 617   numpts = 0;
 618   s[0] = 0;
 619   t[0] = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 620   numpts++;
 621   while (iter) {
 622     s[numpts] = s[numpts-1] + vik_coord_diff ( &(VIK_TRACKPOINT(iter->prev->data)->coord), &(VIK_TRACKPOINT(iter->data)->coord) );
 623     t[numpts] = VIK_TRACKPOINT(iter->data)->timestamp;
 624     numpts++;
 625     iter = iter->next;
 626   }
 627
 628   /* In the following computation, we iterate through periods of time of duration chunk_dur.
 629    * The first period begins at the beginning of the track.  The last period ends at the end of the track.
 630    */
 631   index = 0; /* index of the current trackpoint. */
 632   for (i = 0; i < num_chunks; i++) {
 633     /* we are now covering the interval from t[0] + i*chunk_dur to t[0] + (i+1)*chunk_dur.
 634      * find the first trackpoint outside the current interval, averaging the speeds between intermediate trackpoints.
 635      */
 636     if (t[0] + i*chunk_dur >= t[index]) {
 637       gdouble acc_t = 0, acc_s = 0;
 638       while (t[0] + i*chunk_dur >= t[index]) {
 639         acc_s += (s[index+1]-s[index]);
 640         acc_t += (t[index+1]-t[index]);
 641         index++;
 642       }
 643       v[i] = acc_s/acc_t;
 644     }
 645     else if (i) {
 646       v[i] = v[i-1];
 647     }
 648     else {
 649       v[i] = 0;
 650     }
 651   }
 652   g_free(s);
 653   g_free(t);
 654   return v;
 655 }
 656
 657 /**
 658  * Make a distance/time map, heavily based on the vik_track_make_speed_map method
 659  */
 660 gdouble *vik_track_make_distance_map ( const VikTrack *tr, guint16 num_chunks )
 661 {
 662   gdouble *v, *s, *t;
 663   gdouble duration, chunk_dur;
 664   time_t t1, t2;
 665   int i, pt_count, numpts, index;
 666   GList *iter;
 667
 668   if ( ! tr->trackpoints )
 669     return NULL;
 670
 671   t1 = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 672   t2 = VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->timestamp;
 673   duration = t2 - t1;
 674
 675   if ( !t1 || !t2 || !duration )
 676     return NULL;
 677
 678   if (duration < 0) {
 679     g_warning("negative duration: unsorted trackpoint timestamps?");
 680     return NULL;
 681   }
 682   pt_count = vik_track_get_tp_count(tr);
 683
 684   v = g_malloc ( sizeof(gdouble) * num_chunks );
 685   chunk_dur = duration / num_chunks;
 686
 687   s = g_malloc(sizeof(double) * pt_count);
 688   t = g_malloc(sizeof(double) * pt_count);
 689
 690   iter = tr->trackpoints->next;
 691   numpts = 0;
 692   s[0] = 0;
 693   t[0] = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 694   numpts++;
 695   while (iter) {
 696     s[numpts] = s[numpts-1] + vik_coord_diff ( &(VIK_TRACKPOINT(iter->prev->data)->coord), &(VIK_TRACKPOINT(iter->data)->coord) );
 697     t[numpts] = VIK_TRACKPOINT(iter->data)->timestamp;
 698     numpts++;
 699     iter = iter->next;
 700   }
 701
 702   /* In the following computation, we iterate through periods of time of duration chunk_dur.
 703    * The first period begins at the beginning of the track.  The last period ends at the end of the track.
 704    */
 705   index = 0; /* index of the current trackpoint. */
 706   for (i = 0; i < num_chunks; i++) {
 707     /* we are now covering the interval from t[0] + i*chunk_dur to t[0] + (i+1)*chunk_dur.
 708      * find the first trackpoint outside the current interval, averaging the distance between intermediate trackpoints.
 709      */
 710     if (t[0] + i*chunk_dur >= t[index]) {
 711       gdouble acc_s = 0; // No need for acc_t
 712       while (t[0] + i*chunk_dur >= t[index]) {
 713         acc_s += (s[index+1]-s[index]);
 714         index++;
 715       }
 716       // The only bit that's really different from the speed map - just keep an accululative record distance
 717       v[i] = i ? v[i-1]+acc_s : acc_s;
 718     }
 719     else if (i) {
 720       v[i] = v[i-1];
 721     }
 722     else {
 723       v[i] = 0;
 724     }
 725   }
 726   g_free(s);
 727   g_free(t);
 728   return v;
 729 }
 730
 731 /**
 732  * This uses the 'time' based method to make the graph, (which is a simpler compared to the elevation/distance)
 733  * This results in a slightly blocky graph when it does not have many trackpoints: <60
 734  * NB Somehow the elevation/distance applies some kind of smoothing algorithm,
 735  *   but I don't think any one understands it any more (I certainly don't ATM)
 736  */
 737 gdouble *vik_track_make_elevation_time_map ( const VikTrack *tr, guint16 num_chunks )
 738 {
 739   time_t t1, t2;
 740   gdouble duration, chunk_dur;
 741   GList *iter = tr->trackpoints;
 742
 743   if (!iter || !iter->next) /* zero- or one-point track */
 744     return NULL;
 745
 746   /* test if there's anything worth calculating */
 747   gboolean okay = FALSE;
 748   while ( iter ) {
 749     if ( VIK_TRACKPOINT(iter->data)->altitude != VIK_DEFAULT_ALTITUDE ) {
 750       okay = TRUE;
 751       break;
 752     }
 753     iter = iter->next;
 754   }
 755   if ( ! okay )
 756     return NULL;
 757
 758   t1 = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 759   t2 = VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->timestamp;
 760   duration = t2 - t1;
 761
 762   if ( !t1 || !t2 || !duration )
 763     return NULL;
 764
 765   if (duration < 0) {
 766     g_warning("negative duration: unsorted trackpoint timestamps?");
 767     return NULL;
 768   }
 769   gint pt_count = vik_track_get_tp_count(tr);
 770
 771   // Reset iterator back to the beginning
 772   iter = tr->trackpoints;
 773
 774   gdouble *pts = g_malloc ( sizeof(gdouble) * num_chunks ); // The return altitude values
 775   gdouble *s = g_malloc(sizeof(double) * pt_count); // calculation altitudes
 776   gdouble *t = g_malloc(sizeof(double) * pt_count); // calculation times
 777
 778   chunk_dur = duration / num_chunks;
 779
 780   s[0] = VIK_TRACKPOINT(iter->data)->altitude;
 781   t[0] = VIK_TRACKPOINT(iter->data)->timestamp;
 782   iter = tr->trackpoints->next;
 783   gint numpts = 1;
 784   while (iter) {
 785     s[numpts] = VIK_TRACKPOINT(iter->data)->altitude;
 786     t[numpts] = VIK_TRACKPOINT(iter->data)->timestamp;
 787     numpts++;
 788     iter = iter->next;
 789   }
 790
 791  /* In the following computation, we iterate through periods of time of duration chunk_dur.
 792    * The first period begins at the beginning of the track.  The last period ends at the end of the track.
 793    */
 794   gint index = 0; /* index of the current trackpoint. */
 795   gint i;
 796   for (i = 0; i < num_chunks; i++) {
 797     /* we are now covering the interval from t[0] + i*chunk_dur to t[0] + (i+1)*chunk_dur.
 798      * find the first trackpoint outside the current interval, averaging the heights between intermediate trackpoints.
 799      */
 800     if (t[0] + i*chunk_dur >= t[index]) {
 801       gdouble acc_s = s[index]; // initialise to first point
 802       while (t[0] + i*chunk_dur >= t[index]) {
 803         acc_s += (s[index+1]-s[index]);
 804         index++;
 805       }
 806       pts[i] = acc_s;
 807     }
 808     else if (i) {
 809       pts[i] = pts[i-1];
 810     }
 811     else {
 812       pts[i] = 0;
 813     }
 814   }
 815   g_free(s);
 816   g_free(t);
 817
 818   return pts;
 819 }
 820
 821 /**
 822  * Make a speed/distance map
 823  */
 824 gdouble *vik_track_make_speed_dist_map ( const VikTrack *tr, guint16 num_chunks )
 825 {
 826   gdouble *v, *s, *t;
 827   time_t t1, t2;
 828   gint i, pt_count, numpts, index;
 829   GList *iter;
 830   gdouble duration, total_length, chunk_length;
 831
 832   if ( ! tr->trackpoints )
 833     return NULL;
 834
 835   t1 = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 836   t2 = VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->timestamp;
 837   duration = t2 - t1;
 838
 839   if ( !t1 || !t2 || !duration )
 840     return NULL;
 841
 842   if (duration < 0) {
 843     g_warning("negative duration: unsorted trackpoint timestamps?");
 844     return NULL;
 845   }
 846
 847   total_length = vik_track_get_length_including_gaps ( tr );
 848   chunk_length = total_length / num_chunks;
 849   pt_count = vik_track_get_tp_count(tr);
 850
 851   if (chunk_length <= 0) {
 852     return NULL;
 853   }
 854
 855   v = g_malloc ( sizeof(gdouble) * num_chunks );
 856   s = g_malloc ( sizeof(double) * pt_count );
 857   t = g_malloc ( sizeof(double) * pt_count );
 858
 859   // No special handling of segments ATM...
 860   iter = tr->trackpoints->next;
 861   numpts = 0;
 862   s[0] = 0;
 863   t[0] = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 864   numpts++;
 865   while (iter) {
 866     s[numpts] = s[numpts-1] + vik_coord_diff ( &(VIK_TRACKPOINT(iter->prev->data)->coord), &(VIK_TRACKPOINT(iter->data)->coord) );
 867     t[numpts] = VIK_TRACKPOINT(iter->data)->timestamp;
 868     numpts++;
 869     iter = iter->next;
 870   }
 871
 872   // Iterate through a portion of the track to get an average speed for that part
 873   // This will essentially interpolate between segments, which I think is right given the usage of 'get_length_including_gaps'
 874   index = 0; /* index of the current trackpoint. */
 875   for (i = 0; i < num_chunks; i++) {
 876     // Similar to the make_speed_map, but instead of using a time chunk, use a distance chunk
 877     if (s[0] + i*chunk_length >= s[index]) {
 878       gdouble acc_t = 0, acc_s = 0;
 879       while (s[0] + i*chunk_length >= s[index]) {
 880         acc_s += (s[index+1]-s[index]);
 881         acc_t += (t[index+1]-t[index]);
 882         index++;
 883       }
 884       v[i] = acc_s/acc_t;
 885     }
 886     else if (i) {
 887       v[i] = v[i-1];
 888     }
 889     else {
 890       v[i] = 0;
 891     }
 892   }
 893   g_free(s);
 894   g_free(t);
 895   return v;
 896 }
 897
 898 /* by Alex Foobarian */
 899 VikTrackpoint *vik_track_get_closest_tp_by_percentage_dist ( VikTrack *tr, gdouble reldist, gdouble *meters_from_start )
 900 {
 901   gdouble dist = vik_track_get_length_including_gaps(tr) * reldist;
 902   gdouble current_dist = 0.0;
 903   gdouble current_inc = 0.0;
 904   if ( tr->trackpoints )
 905   {
 906     GList *iter = tr->trackpoints->next;
 907     GList *last_iter = NULL;
 908     gdouble last_dist = 0.0;
 909     while (iter)
 910     {
 911       current_inc = vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
 912                                      &(VIK_TRACKPOINT(iter->prev->data)->coord) );
 913       last_dist = current_dist;
 914       current_dist += current_inc;
 915       if ( current_dist >= dist )
 916         break;
 917       last_iter = iter;
 918       iter = iter->next;
 919     }
 920     if (!iter) { /* passing the end the track */
 921       if (last_iter) {
 922         if (meters_from_start)
 923           *meters_from_start = last_dist;
 924         return(VIK_TRACKPOINT(last_iter->data));
 925       }
 926       else
 927         return NULL;
 928     }
 929     /* we've gone past the dist already, was prev trackpoint closer? */
 930     /* should do a vik_coord_average_weighted() thingy. */
 931     if ( iter->prev && abs(current_dist-current_inc-dist) < abs(current_dist-dist) ) {
 932       if (meters_from_start)
 933         *meters_from_start = last_dist;
 934       iter = iter->prev;
 935     }
 936     else
 937       if (meters_from_start)
 938         *meters_from_start = current_dist;
 939
 940     return VIK_TRACKPOINT(iter->data);
 941
 942   }
 943   return NULL;
 944 }
 945
 946 VikTrackpoint *vik_track_get_closest_tp_by_percentage_time ( VikTrack *tr, gdouble reltime, time_t *seconds_from_start )
 947 {
 948   time_t t_pos, t_start, t_end, t_total;
 949   t_start = VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 950   t_end = VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->timestamp;
 951   t_total = t_end - t_start;
 952
 953   t_pos = t_start + t_total * reltime;
 954
 955   if ( !tr->trackpoints )
 956     return NULL;
 957
 958   GList *iter = tr->trackpoints;
 959
 960   while (iter) {
 961     if (VIK_TRACKPOINT(iter->data)->timestamp == t_pos)
 962       break;
 963     if (VIK_TRACKPOINT(iter->data)->timestamp > t_pos) {
 964       if (iter->prev == NULL)  /* first trackpoint */
 965         break;
 966       time_t t_before = t_pos - VIK_TRACKPOINT(iter->prev)->timestamp;
 967       time_t t_after = VIK_TRACKPOINT(iter->data)->timestamp - t_pos;
 968       if (t_before <= t_after)
 969         iter = iter->prev;
 970       break;
 971     }
 972     else if ((iter->next == NULL) && (t_pos < (VIK_TRACKPOINT(iter->data)->timestamp + 3))) /* last trackpoint: accommodate for round-off */
 973       break;
 974     iter = iter->next;
 975   }
 976
 977   if (!iter)
 978     return NULL;
 979   if (seconds_from_start)
 980     *seconds_from_start = VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(tr->trackpoints->data)->timestamp;
 981   return VIK_TRACKPOINT(iter->data);
 982 }
 983
 984 VikTrackpoint* vik_track_get_tp_by_max_speed ( const VikTrack *tr )
 985 {
 986   gdouble maxspeed = 0.0, speed = 0.0;
 987
 988   if ( !tr->trackpoints )
 989     return NULL;
 990
 991   GList *iter = tr->trackpoints;
 992   VikTrackpoint *max_speed_tp = NULL;
 993
 994   while (iter) {
 995     if (iter->prev) {
 996       if ( VIK_TRACKPOINT(iter->data)->has_timestamp &&
 997            VIK_TRACKPOINT(iter->prev->data)->has_timestamp &&
 998            (! VIK_TRACKPOINT(iter->data)->newsegment) ) {
 999         speed =  vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord), &(VIK_TRACKPOINT(iter->prev->data)->coord) )
1000           / ABS(VIK_TRACKPOINT(iter->data)->timestamp - VIK_TRACKPOINT(iter->prev->data)->timestamp);
1001         if ( speed > maxspeed ) {
1002           maxspeed = speed;
1003           max_speed_tp = VIK_TRACKPOINT(iter->data);
1004         }
1005       }
1006     }
1007     iter = iter->next;
1008   }
1009
1010   if (!max_speed_tp)
1011     return NULL;
1012
1013   return max_speed_tp;
1014 }
1015
1016 VikTrackpoint* vik_track_get_tp_by_max_alt ( const VikTrack *tr )
1017 {
1018   gdouble maxalt = -5000.0;
1019   if ( !tr->trackpoints )
1020     return NULL;
1021
1022   GList *iter = tr->trackpoints;
1023   VikTrackpoint *max_alt_tp = NULL;
1024
1025   while (iter) {
1026     if ( VIK_TRACKPOINT(iter->data)->altitude > maxalt ) {
1027       maxalt = VIK_TRACKPOINT(iter->data)->altitude;
1028       max_alt_tp = VIK_TRACKPOINT(iter->data);
1029     }
1030     iter = iter->next;
1031   }
1032
1033   if (!max_alt_tp)
1034     return NULL;
1035
1036   return max_alt_tp;
1037 }
1038
1039 VikTrackpoint* vik_track_get_tp_by_min_alt ( const VikTrack *tr )
1040 {
1041   gdouble minalt = 25000.0;
1042   if ( !tr->trackpoints )
1043     return NULL;
1044
1045   GList *iter = tr->trackpoints;
1046   VikTrackpoint *min_alt_tp = NULL;
1047
1048   while (iter) {
1049     if ( VIK_TRACKPOINT(iter->data)->altitude < minalt ) {
1050       minalt = VIK_TRACKPOINT(iter->data)->altitude;
1051       min_alt_tp = VIK_TRACKPOINT(iter->data);
1052     }
1053     iter = iter->next;
1054   }
1055
1056   if (!min_alt_tp)
1057     return NULL;
1058
1059   return min_alt_tp;
1060 }
1061
1062 gboolean vik_track_get_minmax_alt ( const VikTrack *tr, gdouble *min_alt, gdouble *max_alt )
1063 {
1064   *min_alt = 25000;
1065   *max_alt = -5000;
1066   if ( tr && tr->trackpoints && tr->trackpoints->data && (VIK_TRACKPOINT(tr->trackpoints->data)->altitude != VIK_DEFAULT_ALTITUDE) ) {
1067     GList *iter = tr->trackpoints->next;
1068     gdouble tmp_alt;
1069     while (iter)
1070     {
1071       tmp_alt = VIK_TRACKPOINT(iter->data)->altitude;
1072       if ( tmp_alt > *max_alt )
1073         *max_alt = tmp_alt;
1074       if ( tmp_alt < *min_alt )
1075         *min_alt = tmp_alt;
1076       iter = iter->next;
1077     }
1078     return TRUE;
1079   }
1080   return FALSE;
1081 }
1082
1083 void vik_track_marshall ( VikTrack *tr, guint8 **data, guint *datalen)
1084 {
1085   GList *tps;
1086   GByteArray *b = g_byte_array_new();
1087   guint len;
1088   guint intp, ntp;
1089
1090   g_byte_array_append(b, (guint8 *)tr, sizeof(*tr));
1091
1092   /* we'll fill out number of trackpoints later */
1093   intp = b->len;
1094   g_byte_array_append(b, (guint8 *)&len, sizeof(len));
1095
1096   tps = tr->trackpoints;
1097   ntp = 0;
1098   while (tps) {
1099     g_byte_array_append(b, (guint8 *)tps->data, sizeof(VikTrackpoint));
1100     tps = tps->next;
1101     ntp++;
1102   }
1103   *(guint *)(b->data + intp) = ntp;
1104
1105   len = (tr->name) ? strlen(tr->name)+1 : 0;
1106   g_byte_array_append(b, (guint8 *)&len, sizeof(len));
1107   if (tr->name) g_byte_array_append(b, (guint8 *)tr->name, len);
1108
1109   len = (tr->comment) ? strlen(tr->comment)+1 : 0;
1110   g_byte_array_append(b, (guint8 *)&len, sizeof(len));
1111   if (tr->comment) g_byte_array_append(b, (guint8 *)tr->comment, len);
1112
1113   *data = b->data;
1114   *datalen = b->len;
1115   g_byte_array_free(b, FALSE);
1116 }
1117
1118 VikTrack *vik_track_unmarshall (guint8 *data, guint datalen)
1119 {
1120   guint len;
1121   VikTrack *new_tr = vik_track_new();
1122   VikTrackpoint *new_tp;
1123   guint ntp;
1124   gint i;
1125
1126   /* only the visibility is needed */
1127   new_tr->visible = ((VikTrack *)data)->visible;
1128   data += sizeof(*new_tr);
1129
1130   ntp = *(guint *)data;
1131   data += sizeof(ntp);
1132
1133   for (i=0; i<ntp; i++) {
1134     new_tp = vik_trackpoint_new();
1135     memcpy(new_tp, data, sizeof(*new_tp));
1136     data += sizeof(*new_tp);
1137     new_tr->trackpoints = g_list_append(new_tr->trackpoints, new_tp);
1138   }
1139
1140   len = *(guint *)data;
1141   data += sizeof(len);
1142   if (len) {
1143     new_tr->name = g_strdup((gchar *)data);
1144   }
1145   data += len;
1146
1147   len = *(guint *)data;
1148   data += sizeof(len);
1149   if (len) {
1150     new_tr->comment = g_strdup((gchar *)data);
1151   }
1152   return new_tr;
1153 }
1154
1155 void vik_track_apply_dem_data ( VikTrack *tr )
1156 {
1157   GList *tp_iter;
1158   gint16 elev;
1159   tp_iter = tr->trackpoints;
1160   while ( tp_iter ) {
1161     /* TODO: of the 4 possible choices we have for choosing an elevation
1162      * (trackpoint in between samples), choose the one with the least elevation change
1163      * as the last */
1164     elev = a_dems_get_elev_by_coord ( &(VIK_TRACKPOINT(tp_iter->data)->coord), VIK_DEM_INTERPOL_BEST );
1165     if ( elev != VIK_DEM_INVALID_ELEVATION )
1166       VIK_TRACKPOINT(tp_iter->data)->altitude = elev;
1167     tp_iter = tp_iter->next;
1168   }
1169 }
1170
1171 /*
1172  * Apply DEM data (if available) - to only the last trackpoint
1173  */
1174 void vik_track_apply_dem_data_last_trackpoint ( VikTrack *tr )
1175 {
1176   gint16 elev;
1177   if ( tr->trackpoints ) {
1178     /* As in vik_track_apply_dem_data above - use 'best' interpolation method */
1179     elev = a_dems_get_elev_by_coord ( &(VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->coord), VIK_DEM_INTERPOL_BEST );
1180     if ( elev != VIK_DEM_INVALID_ELEVATION )
1181       VIK_TRACKPOINT(g_list_last(tr->trackpoints)->data)->altitude = elev;
1182   }
1183 }
1184
1185 /* appends t2 to t1, leaving t2 with no trackpoints */
1186 void vik_track_steal_and_append_trackpoints ( VikTrack *t1, VikTrack *t2 )
1187 {
1188   if ( t1->trackpoints ) {
1189     GList *tpiter = t1->trackpoints;
1190     while ( tpiter->next )
1191       tpiter = tpiter->next;
1192     tpiter->next = t2->trackpoints;
1193     t2->trackpoints->prev = tpiter;
1194   } else
1195     t1->trackpoints = t2->trackpoints;
1196   t2->trackpoints = NULL;
1197 }
1198
1199 /* starting at the end, looks backwards for the last "double point", a duplicate trackpoint.
1200  * If there is no double point, deletes all the trackpoints.
1201  * Returns the new end of the track (or the start if there are no double points)
1202  */
1203 VikCoord *vik_track_cut_back_to_double_point ( VikTrack *tr )
1204 {
1205   GList *iter = tr->trackpoints;
1206   VikCoord *rv;
1207
1208   if ( !iter )
1209     return NULL;
1210   while ( iter->next )
1211     iter = iter->next;
1212
1213
1214   while ( iter->prev ) {
1215     if ( vik_coord_equals((VikCoord *)iter->data, (VikCoord *)iter->prev->data) ) {
1216       GList *prev = iter->prev;
1217
1218       rv = g_malloc(sizeof(VikCoord));
1219       *rv = *((VikCoord *) iter->data);
1220
1221       /* truncate trackpoint list */
1222       iter->prev = NULL; /* pretend it's the end */
1223       g_list_foreach ( iter, (GFunc) g_free, NULL );
1224       g_list_free( iter );
1225
1226       prev->next = NULL;
1227
1228       return rv;
1229     }
1230     iter = iter->prev;
1231   }
1232
1233   /* no double point found! */
1234   rv = g_malloc(sizeof(VikCoord));
1235   *rv = *((VikCoord *) tr->trackpoints->data);
1236   g_list_foreach ( tr->trackpoints, (GFunc) g_free, NULL );
1237   g_list_free( tr->trackpoints );
1238   tr->trackpoints = NULL;
1239   return rv;
1240 }
1241