Routino SVN Repository Browser

Check out the latest version of Routino: svn co http://routino.org/svn/trunk routino

ViewVC logotype

Contents of /trunk/src/prunex.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1999 - (show annotations) (download) (as text)
Sat Jul 27 10:33:04 2019 UTC (5 years, 7 months ago) by amb
File MIME type: text/x-csrc
File size: 42338 byte(s)
Add more checking of memory allocation success/failure by combining
the allocation and the assert into one function.

1 /***************************************
2 Data pruning functions.
3
4 Part of the Routino routing software.
5 ******************/ /******************
6 This file Copyright 2011-2014, 2019 Andrew M. Bishop
7
8 This program is free software: you can redistribute it and/or modify
9 it under the terms of the GNU Affero General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Affero General Public License for more details.
17
18 You should have received a copy of the GNU Affero General Public License
19 along with this program. If not, see <http://www.gnu.org/licenses/>.
20 ***************************************/
21
22
23 #include <stdlib.h>
24
25 #include "types.h"
26 #include "segments.h"
27
28 #include "typesx.h"
29 #include "nodesx.h"
30 #include "segmentsx.h"
31 #include "waysx.h"
32
33 #include "prunex.h"
34
35 #include "files.h"
36 #include "logging.h"
37
38
39 /* Global variables */
40
41 /*+ The command line '--tmpdir' option or its default value. +*/
42 extern char *option_tmpdirname;
43
44 /* Local functions */
45
46 static void prune_segment(SegmentsX *segmentsx,SegmentX *segmentx);
47 static void modify_segment(SegmentsX *segmentsx,SegmentX *segmentx,index_t newnode1,index_t newnode2);
48
49 static void unlink_segment_node1_refs(SegmentsX *segmentsx,SegmentX *segmentx);
50 static void unlink_segment_node2_refs(SegmentsX *segmentsx,SegmentX *segmentx);
51
52 static double distance(double lat1,double lon1,double lat2,double lon2);
53
54
55 /*++++++++++++++++++++++++++++++++++++++
56 Initialise the data structures needed for pruning.
57
58 NodesX *nodesx The set of nodes to use.
59
60 SegmentsX *segmentsx The set of segments to use.
61
62 WaysX *waysx The set of ways to use.
63 ++++++++++++++++++++++++++++++++++++++*/
64
65 void StartPruning(NodesX *nodesx,SegmentsX *segmentsx,WaysX *waysx)
66 {
67 SegmentX segmentx;
68 index_t index=0,lastnode1=NO_NODE;
69
70 if(segmentsx->number==0)
71 return;
72
73 /* Print the start message */
74
75 printf_first("Adding Extra Segment Indexes: Segments=0");
76
77 /* Allocate the array of next segment */
78
79 segmentsx->next1=(index_t*)calloc_logassert(segmentsx->number,sizeof(index_t));
80 log_malloc(segmentsx->next1,segmentsx->number*sizeof(index_t));
81
82 /* Open the file read-only */
83
84 segmentsx->fd=ReOpenFileBuffered(segmentsx->filename_tmp);
85
86 /* Read the on-disk image */
87
88 while(!ReadFileBuffered(segmentsx->fd,&segmentx,sizeof(SegmentX)))
89 {
90 index_t node1=segmentx.node1;
91
92 if(index==0)
93 ;
94 else if(lastnode1==node1)
95 segmentsx->next1[index-1]=index;
96 else
97 segmentsx->next1[index-1]=NO_SEGMENT;
98
99 lastnode1=node1;
100 index++;
101
102 if(!(index%10000))
103 printf_middle("Added Extra Segment Indexes: Segments=%"Pindex_t,index);
104 }
105
106 segmentsx->next1[index-1]=NO_SEGMENT;
107
108 /* Close the file */
109
110 segmentsx->fd=CloseFileBuffered(segmentsx->fd);
111
112 /* Print the final message */
113
114 printf_last("Added Extra Segment Indexes: Segments=%"Pindex_t,segmentsx->number);
115 }
116
117
118 /*++++++++++++++++++++++++++++++++++++++
119 Delete the data structures needed for pruning.
120
121 NodesX *nodesx The set of nodes to use.
122
123 SegmentsX *segmentsx The set of segments to use.
124
125 WaysX *waysx The set of ways to use.
126 ++++++++++++++++++++++++++++++++++++++*/
127
128 void FinishPruning(NodesX *nodesx,SegmentsX *segmentsx,WaysX *waysx)
129 {
130 if(segmentsx->next1)
131 {
132 log_free(segmentsx->next1);
133 free(segmentsx->next1);
134 segmentsx->next1=NULL;
135 }
136 }
137
138
139 /*++++++++++++++++++++++++++++++++++++++
140 Prune out any groups of nodes and segments whose total length is less than a
141 specified minimum.
142
143 NodesX *nodesx The set of nodes to use.
144
145 SegmentsX *segmentsx The set of segments to use.
146
147 WaysX *waysx The set of ways to use.
148
149 distance_t minimum The minimum distance to keep.
150 ++++++++++++++++++++++++++++++++++++++*/
151
152 void PruneIsolatedRegions(NodesX *nodesx,SegmentsX *segmentsx,WaysX *waysx,distance_t minimum)
153 {
154 WaysX *newwaysx;
155 WayX tmpwayx;
156 transport_t transport;
157 BitMask *connected,*region;
158 index_t *regionsegments,*othersegments;
159 index_t nallocregionsegments,nallocothersegments;
160
161 if(nodesx->number==0 || segmentsx->number==0)
162 return;
163
164 /* Map into memory / open the files */
165
166 #if !SLIM
167 nodesx->data=MapFile(nodesx->filename_tmp);
168 segmentsx->data=MapFileWriteable(segmentsx->filename_tmp);
169 waysx->data=MapFile(waysx->filename_tmp);
170 #else
171 nodesx->fd=SlimMapFile(nodesx->filename_tmp);
172 segmentsx->fd=SlimMapFileWriteable(segmentsx->filename_tmp);
173 waysx->fd=SlimMapFile(waysx->filename_tmp);
174
175 InvalidateNodeXCache(nodesx->cache);
176 InvalidateSegmentXCache(segmentsx->cache);
177 InvalidateWayXCache(waysx->cache);
178 #endif
179
180 newwaysx=NewWayList(0,0);
181 CloseFileBuffered(newwaysx->fd);
182
183 newwaysx->fd=SlimMapFileWriteable(newwaysx->filename_tmp);
184
185 connected=AllocBitMask(segmentsx->number);
186 region =AllocBitMask(segmentsx->number);
187
188 log_malloc(connected,LengthBitMask(segmentsx->number)*sizeof(BitMask));
189 log_malloc(region ,LengthBitMask(segmentsx->number)*sizeof(BitMask));
190
191 regionsegments=(index_t*)malloc_logassert((nallocregionsegments=1024)*sizeof(index_t));
192 othersegments =(index_t*)malloc_logassert((nallocothersegments =1024)*sizeof(index_t));
193
194 /* Loop through the transport types */
195
196 for(transport=Transport_None+1;transport<Transport_Count;transport++)
197 {
198 index_t i,j;
199 index_t nregions=0,npruned=0,nadjusted=0;
200 const char *transport_str=TransportName(transport);
201 transports_t transports=TRANSPORTS(transport);
202
203 if(!(waysx->transports&transports))
204 continue;
205
206 /* Print the start message */
207
208 printf_first("Pruning Isolated Regions (%s): Segments=0 Adjusted=0 Pruned=0",transport_str);
209
210 /* Loop through the segments and find the disconnected ones */
211
212 ClearAllBits(connected,segmentsx->number);
213 ClearAllBits(region ,segmentsx->number);
214
215 for(i=0;i<segmentsx->number;i++)
216 {
217 index_t nregionsegments=0,nothersegments=0;
218 distance_t total=0;
219 SegmentX *segmentx;
220 WayX *wayx;
221
222 if(IsBitSet(connected,i))
223 goto endloop;
224
225 segmentx=LookupSegmentX(segmentsx,i,1);
226
227 if(IsPrunedSegmentX(segmentx))
228 goto endloop;
229
230 if(segmentx->way<waysx->number)
231 wayx=LookupWayX(waysx,segmentx->way,1);
232 else
233 SlimFetch(newwaysx->fd,(wayx=&tmpwayx),sizeof(WayX),(segmentx->way-waysx->number)*sizeof(WayX));
234
235 if(!(wayx->way.allow&transports))
236 goto endloop;
237
238 othersegments[nothersegments++]=i;
239 SetBit(region,i);
240
241 do
242 {
243 index_t thissegment,nodes[2];
244
245 thissegment=othersegments[--nothersegments];
246
247 if(nregionsegments==nallocregionsegments)
248 regionsegments=(index_t*)realloc_logassert(regionsegments,(nallocregionsegments+=1024)*sizeof(index_t));
249
250 regionsegments[nregionsegments++]=thissegment;
251
252 segmentx=LookupSegmentX(segmentsx,thissegment,1);
253
254 nodes[0]=segmentx->node1;
255 nodes[1]=segmentx->node2;
256 total+=DISTANCE(segmentx->distance);
257
258 for(j=0;j<2;j++)
259 {
260 NodeX *nodex=LookupNodeX(nodesx,nodes[j],1);
261
262 if(!(nodex->allow&transports))
263 continue;
264
265 segmentx=FirstSegmentX(segmentsx,nodes[j],1);
266
267 while(segmentx)
268 {
269 index_t segment=IndexSegmentX(segmentsx,segmentx);
270
271 if(segment!=thissegment)
272 {
273 if(segmentx->way<waysx->number)
274 wayx=LookupWayX(waysx,segmentx->way,1);
275 else
276 SlimFetch(newwaysx->fd,(wayx=&tmpwayx),sizeof(WayX),(segmentx->way-waysx->number)*sizeof(WayX));
277
278 if(wayx->way.allow&transports)
279 {
280 /* Already connected - finish */
281
282 if(IsBitSet(connected,segment))
283 {
284 total=minimum;
285 goto foundconnection;
286 }
287
288 /* Not in region - add to list */
289
290 if(!IsBitSet(region,segment))
291 {
292 if(nothersegments==nallocothersegments)
293 othersegments=(index_t*)realloc_logassert(othersegments,(nallocothersegments+=1024)*sizeof(index_t));
294
295 othersegments[nothersegments++]=segment;
296 SetBit(region,segment);
297 }
298 }
299 }
300
301 segmentx=NextSegmentX(segmentsx,segmentx,nodes[j]);
302 }
303 }
304 }
305 while(nothersegments>0 && total<minimum);
306
307 foundconnection:
308
309 /* Prune the segments or mark them as connected */
310
311 if(total<minimum) /* not connected - delete them */
312 {
313 nregions++;
314
315 for(j=0;j<nregionsegments;j++)
316 {
317 SegmentX *segmentx;
318 WayX *wayx,tmpwayx;
319
320 SetBit(connected,regionsegments[j]); /* not really connected, but don't need to check again */
321 ClearBit(region,regionsegments[j]);
322
323 segmentx=LookupSegmentX(segmentsx,regionsegments[j],1);
324
325 if(segmentx->way<waysx->number)
326 wayx=LookupWayX(waysx,segmentx->way,1);
327 else
328 SlimFetch(newwaysx->fd,(wayx=&tmpwayx),sizeof(WayX),(segmentx->way-waysx->number)*sizeof(WayX));
329
330 if(wayx->way.allow==transports)
331 {
332 prune_segment(segmentsx,segmentx);
333
334 npruned++;
335 }
336 else
337 {
338 if(segmentx->way<waysx->number) /* create a new way */
339 {
340 tmpwayx=*wayx;
341
342 tmpwayx.way.allow&=~transports;
343
344 segmentx->way=waysx->number+newwaysx->number;
345
346 SlimReplace(newwaysx->fd,&tmpwayx,sizeof(WayX),(segmentx->way-waysx->number)*sizeof(WayX));
347
348 newwaysx->number++;
349
350 PutBackSegmentX(segmentsx,segmentx);
351 }
352 else /* modify the existing one */
353 {
354 tmpwayx.way.allow&=~transports;
355
356 SlimReplace(newwaysx->fd,&tmpwayx,sizeof(WayX),(segmentx->way-waysx->number)*sizeof(WayX));
357 }
358
359 nadjusted++;
360 }
361 }
362 }
363 else /* connected - mark as part of the main region */
364 {
365 for(j=0;j<nregionsegments;j++)
366 {
367 SetBit(connected,regionsegments[j]);
368 ClearBit(region,regionsegments[j]);
369 }
370
371 for(j=0;j<nothersegments;j++)
372 {
373 SetBit(connected,othersegments[j]);
374 ClearBit(region,othersegments[j]);
375 }
376 }
377
378 endloop:
379
380 if(!((i+1)%10000))
381 printf_middle("Pruning Isolated Regions (%s): Segments=%"Pindex_t" Adjusted=%"Pindex_t" Pruned=%"Pindex_t" (%"Pindex_t" Regions)",transport_str,i+1,nadjusted,npruned,nregions);
382 }
383
384 /* Print the final message */
385
386 printf_last("Pruned Isolated Regions (%s): Segments=%"Pindex_t" Adjusted=%"Pindex_t" Pruned=%"Pindex_t" (%"Pindex_t" Regions)",transport_str,segmentsx->number,nadjusted,npruned,nregions);
387 }
388
389 /* Unmap from memory / close the files */
390
391 log_free(region);
392 log_free(connected);
393
394 free(region);
395 free(connected);
396
397 free(regionsegments);
398 free(othersegments);
399
400 #if !SLIM
401 nodesx->data=UnmapFile(nodesx->data);
402 segmentsx->data=UnmapFile(segmentsx->data);
403 waysx->data=UnmapFile(waysx->data);
404 #else
405 nodesx->fd=SlimUnmapFile(nodesx->fd);
406 segmentsx->fd=SlimUnmapFile(segmentsx->fd);
407 waysx->fd=SlimUnmapFile(waysx->fd);
408 #endif
409
410 SlimUnmapFile(newwaysx->fd);
411
412 waysx->number+=newwaysx->number;
413
414 waysx->fd=OpenFileBufferedAppend(waysx->filename_tmp);
415
416 newwaysx->fd=ReOpenFileBuffered(newwaysx->filename_tmp);
417
418 while(!ReadFileBuffered(newwaysx->fd,&tmpwayx,sizeof(WayX)))
419 WriteFileBuffered(waysx->fd,&tmpwayx,sizeof(WayX));
420
421 CloseFileBuffered(waysx->fd);
422 CloseFileBuffered(newwaysx->fd);
423
424 FreeWayList(newwaysx,0);
425 }
426
427
428 /*++++++++++++++++++++++++++++++++++++++
429 Prune out any segments that are shorter than a specified minimum.
430
431 NodesX *nodesx The set of nodes to use.
432
433 SegmentsX *segmentsx The set of segments to use.
434
435 WaysX *waysx The set of ways to use.
436
437 distance_t minimum The maximum length to remove or one less than the minimum length to keep.
438 ++++++++++++++++++++++++++++++++++++++*/
439
440 void PruneShortSegments(NodesX *nodesx,SegmentsX *segmentsx,WaysX *waysx,distance_t minimum)
441 {
442 index_t i;
443 index_t nshort=0,npruned=0;
444
445 if(nodesx->number==0 || segmentsx->number==0 || waysx->number==0)
446 return;
447
448 /* Print the start message */
449
450 printf_first("Pruning Short Segments: Segments=0 Short=0 Pruned=0");
451
452 /* Map into memory / open the files */
453
454 #if !SLIM
455 nodesx->data=MapFileWriteable(nodesx->filename_tmp);
456 segmentsx->data=MapFileWriteable(segmentsx->filename_tmp);
457 waysx->data=MapFile(waysx->filename_tmp);
458 #else
459 nodesx->fd=SlimMapFileWriteable(nodesx->filename_tmp);
460 segmentsx->fd=SlimMapFileWriteable(segmentsx->filename_tmp);
461 waysx->fd=SlimMapFile(waysx->filename_tmp);
462
463 InvalidateNodeXCache(nodesx->cache);
464 InvalidateSegmentXCache(segmentsx->cache);
465 InvalidateWayXCache(waysx->cache);
466 #endif
467
468 /* Loop through the segments and find the short ones for possible modification */
469
470 for(i=0;i<segmentsx->number;i++)
471 {
472 SegmentX *segmentx2=LookupSegmentX(segmentsx,i,2);
473
474 if(IsPrunedSegmentX(segmentx2))
475 goto endloop;
476
477 /*
478 :
479 Initial state: ..N3 -------- N2
480 : S2
481
482 :
483 Final state: ..N3
484 :
485
486 = OR =
487
488 : :
489 Initial state: ..N1 -------- N2 ---- N3 -------- N4..
490 : S1 S2 S3 :
491
492 : :
493 Final state: ..N1 ------------ N3 ------------ N4..
494 : S1 S3 :
495
496 Not if N1 is the same as N4.
497 Must not delete N2 (or N3) if S2 (or S3) has different one-way properties from S1.
498 Must not delete N2 (or N3) if S2 (or S3) has different highway properties from S1.
499 Must combine N2, S2 and N3 disallowed transports into new N3.
500 Must not delete N2 (or N3) if it is a mini-roundabout.
501 Must not delete N2 (or N3) if it is involved in a turn restriction.
502
503 = OR =
504
505 : :
506 Initial state: ..N1 -------- N2 ---- N3..
507 : S1 S2 :
508
509 : :
510 Final state: ..N1 ------------ N3..
511 : S1 :
512
513 Not if N1 is the same as N3.
514 Not if S1 has different one-way properties from S2.
515 Not if S1 has different highway properties from S2.
516 Not if N2 disallows transports allowed on S1 and S2.
517 Not if N2 is a mini-roundabout.
518 Not if N2 is involved in a turn restriction.
519 */
520
521 if(DISTANCE(segmentx2->distance)<=minimum)
522 {
523 index_t node1=NO_NODE,node2,node3,node4=NO_NODE;
524 index_t segment1=NO_SEGMENT,segment2=i,segment3=NO_SEGMENT;
525 SegmentX *segmentx;
526 int segcount2=0,segcount3=0;
527
528 nshort++;
529
530 node2=segmentx2->node1;
531 node3=segmentx2->node2;
532
533 /* Count the segments connected to N2 */
534
535 segmentx=FirstSegmentX(segmentsx,node2,4);
536
537 while(segmentx)
538 {
539 segcount2++;
540
541 if(segment1==NO_SEGMENT)
542 {
543 index_t segment=IndexSegmentX(segmentsx,segmentx);
544
545 if(segment!=segment2)
546 {
547 segment1=segment;
548 node1=OtherNode(segmentx,node2);
549 }
550 }
551 else if(segcount2>2)
552 break;
553
554 segmentx=NextSegmentX(segmentsx,segmentx,node2);
555 }
556
557 /* Count the segments connected to N3 */
558
559 segmentx=FirstSegmentX(segmentsx,node3,4);
560
561 while(segmentx)
562 {
563 segcount3++;
564
565 if(segment3==NO_SEGMENT)
566 {
567 index_t segment=IndexSegmentX(segmentsx,segmentx);
568
569 if(segment!=segment2)
570 {
571 segment3=segment;
572 node4=OtherNode(segmentx,node3);
573 }
574 }
575 else if(segcount3>2)
576 break;
577
578 segmentx=NextSegmentX(segmentsx,segmentx,node3);
579 }
580
581 /* Check which case we are handling (and canonicalise) */
582
583 if(segcount2>2 && segcount3>2) /* none of the cases in diagram - too complicated */
584 {
585 goto endloop;
586 }
587 else if(segcount2==1 || segcount3==1) /* first case in diagram - prune segment */
588 {
589 prune_segment(segmentsx,segmentx2);
590 }
591 else if(segcount2==2 && segcount3==2) /* second case in diagram - modify one segment and prune segment */
592 {
593 SegmentX *segmentx1,*segmentx3;
594 WayX *wayx1,*wayx2,*wayx3;
595 NodeX *nodex2,*nodex3,*newnodex;
596 index_t newnode;
597 int join12=1,join23=1,same13=1;
598
599 /* Check if pruning would collapse a loop */
600
601 if(node1==node4)
602 goto endloop;
603
604 /* Check if allowed due to one-way properties */
605
606 segmentx1=LookupSegmentX(segmentsx,segment1,1);
607 segmentx3=LookupSegmentX(segmentsx,segment3,3);
608
609 if(!IsOneway(segmentx1) && !IsOneway(segmentx2))
610 ;
611 else if(IsOneway(segmentx1) && IsOneway(segmentx2))
612 {
613 if(IsOnewayTo(segmentx1,node2) && !IsOnewayFrom(segmentx2,node2)) /* S1 is one-way but S2 doesn't continue */
614 join12=0;
615
616 if(IsOnewayFrom(segmentx1,node2) && !IsOnewayTo(segmentx2,node2)) /* S1 is one-way but S2 doesn't continue */
617 join12=0;
618 }
619 else
620 join12=0;
621
622 if(!IsOneway(segmentx3) && !IsOneway(segmentx2))
623 ;
624 else if(IsOneway(segmentx3) && IsOneway(segmentx2))
625 {
626 if(IsOnewayTo(segmentx3,node3) && !IsOnewayFrom(segmentx2,node3)) /* S3 is one-way but S2 doesn't continue */
627 join23=0;
628
629 if(IsOnewayFrom(segmentx3,node3) && !IsOnewayTo(segmentx2,node3)) /* S3 is one-way but S2 doesn't continue */
630 join23=0;
631 }
632 else
633 join23=0;
634
635 if(!join12 && !join23)
636 goto endloop;
637
638 /* Check if allowed due to highway properties */
639
640 wayx1=LookupWayX(waysx,segmentx1->way,1);
641 wayx2=LookupWayX(waysx,segmentx2->way,2);
642 wayx3=LookupWayX(waysx,segmentx3->way,3);
643
644 if(WaysCompare(&wayx1->way,&wayx2->way))
645 join12=0;
646
647 if(WaysCompare(&wayx3->way,&wayx2->way))
648 join23=0;
649
650 if(!join12 && !join23)
651 goto endloop;
652
653 /* Check if allowed due to mini-roundabout and turn restriction */
654
655 nodex2=LookupNodeX(nodesx,node2,2);
656 nodex3=LookupNodeX(nodesx,node3,3);
657
658 if(nodex2->flags&NODE_MINIRNDBT)
659 join12=0;
660
661 if(nodex3->flags&NODE_MINIRNDBT)
662 join23=0;
663
664 if(!join12 && !join23)
665 goto endloop;
666
667 if(nodex2->flags&NODE_TURNRSTRCT2 || nodex2->flags&NODE_TURNRSTRCT)
668 join12=0;
669
670 if(nodex3->flags&NODE_TURNRSTRCT2 || nodex3->flags&NODE_TURNRSTRCT)
671 join23=0;
672
673 if(!join12 && !join23)
674 goto endloop;
675
676 /* New node properties */
677
678 if(join12)
679 {
680 newnode=node3;
681 newnodex=nodex3;
682 }
683 else /* if(join23) */
684 {
685 newnode=node2;
686 newnodex=nodex2;
687 }
688
689 newnodex->allow=nodex2->allow&nodex3->allow; /* combine the restrictions of the two nodes */
690 newnodex->allow&=~((~wayx2->way.allow)&wayx3->way.allow); /* disallow anything blocked by segment2 */
691 newnodex->allow&=~((~wayx2->way.allow)&wayx1->way.allow); /* disallow anything blocked by segment2 */
692
693 newnodex->latitude =(nodex2->latitude +nodex3->latitude )/2;
694 newnodex->longitude=(nodex2->longitude+nodex3->longitude)/2;
695
696 PutBackNodeX(nodesx,newnodex);
697
698 /* Modify segments - update the distances */
699
700 if(!IsOneway(segmentx1) && !IsOneway(segmentx3))
701 ;
702 else if(IsOneway(segmentx1) && IsOneway(segmentx3))
703 {
704 if(IsOnewayTo(segmentx1,node3) && !IsOnewayFrom(segmentx3,node3)) /* S1 is one-way but S3 doesn't continue */
705 same13=0;
706
707 if(IsOnewayFrom(segmentx1,node3) && !IsOnewayTo(segmentx3,node3)) /* S1 is one-way but S3 doesn't continue */
708 same13=0;
709 }
710 else
711 same13=0;
712
713 if(WaysCompare(&wayx1->way,&wayx3->way))
714 same13=0;
715
716 if(same13)
717 {
718 segmentx1->distance+=DISTANCE(segmentx2->distance)/2;
719 segmentx3->distance+=DISTANCE(segmentx2->distance)-DISTANCE(segmentx2->distance)/2;
720 }
721 else if(join12)
722 segmentx1->distance+=DISTANCE(segmentx2->distance);
723 else /* if(join23) */
724 segmentx3->distance+=DISTANCE(segmentx2->distance);
725
726 /* Modify segments - update the segments */
727
728 if(segmentx1->node1==node1)
729 {
730 if(segmentx1->node2!=newnode)
731 modify_segment(segmentsx,segmentx1,node1,newnode);
732 else
733 PutBackSegmentX(segmentsx,segmentx1);
734 }
735 else /* if(segmentx1->node2==node1) */
736 {
737 if(segmentx1->node1!=newnode)
738 modify_segment(segmentsx,segmentx1,newnode,node1);
739 else
740 PutBackSegmentX(segmentsx,segmentx1);
741 }
742
743 if(segmentx3->node1==node4)
744 {
745 if(segmentx3->node2!=newnode)
746 modify_segment(segmentsx,segmentx3,node4,newnode);
747 else
748 PutBackSegmentX(segmentsx,segmentx3);
749 }
750 else /* if(segmentx3->node2==node4) */
751 {
752 if(segmentx3->node1!=newnode)
753 modify_segment(segmentsx,segmentx3,newnode,node4);
754 else
755 PutBackSegmentX(segmentsx,segmentx3);
756 }
757
758 ReLookupSegmentX(segmentsx,segmentx2);
759
760 prune_segment(segmentsx,segmentx2);
761 }
762 else /* third case in diagram - prune one segment */
763 {
764 SegmentX *segmentx1;
765 WayX *wayx1,*wayx2;
766 NodeX *nodex2;
767
768 if(segcount3==2) /* not as in diagram, shuffle things round */
769 {
770 index_t temp;
771
772 temp=segment1; segment1=segment3; segment3=temp;
773 temp=node1; node1=node4; node4=temp;
774 temp=node2; node2=node3; node3=temp;
775 }
776
777 /* Check if pruning would collapse a loop */
778
779 if(node1==node3)
780 goto endloop;
781
782 /* Check if allowed due to one-way properties */
783
784 segmentx1=LookupSegmentX(segmentsx,segment1,1);
785
786 if(!IsOneway(segmentx1) && !IsOneway(segmentx2))
787 ;
788 else if(IsOneway(segmentx1) && IsOneway(segmentx2))
789 {
790 if(IsOnewayTo(segmentx1,node2) && !IsOnewayFrom(segmentx2,node2)) /* S1 is one-way but S2 doesn't continue */
791 goto endloop;
792
793 if(IsOnewayFrom(segmentx1,node2) && !IsOnewayTo(segmentx2,node2)) /* S1 is one-way but S2 doesn't continue */
794 goto endloop;
795 }
796 else
797 goto endloop;
798
799 /* Check if allowed due to highway properties */
800
801 wayx1=LookupWayX(waysx,segmentx1->way,1);
802 wayx2=LookupWayX(waysx,segmentx2->way,2);
803
804 if(WaysCompare(&wayx1->way,&wayx2->way))
805 goto endloop;
806
807 /* Check if allowed due to mini-roundabout and turn restriction */
808
809 nodex2=LookupNodeX(nodesx,node2,2);
810
811 if(nodex2->flags&NODE_MINIRNDBT)
812 goto endloop;
813
814 if(nodex2->flags&NODE_TURNRSTRCT2 || nodex2->flags&NODE_TURNRSTRCT)
815 goto endloop;
816
817 /* Check if allowed due to node restrictions */
818
819 if((nodex2->allow&wayx1->way.allow)!=wayx1->way.allow)
820 goto endloop;
821
822 if((nodex2->allow&wayx2->way.allow)!=wayx2->way.allow)
823 goto endloop;
824
825 /* Modify segments */
826
827 segmentx1->distance+=DISTANCE(segmentx2->distance);
828
829 if(segmentx1->node1==node1)
830 modify_segment(segmentsx,segmentx1,node1,node3);
831 else /* if(segmentx1->node2==node1) */
832 modify_segment(segmentsx,segmentx1,node3,node1);
833
834 ReLookupSegmentX(segmentsx,segmentx2);
835
836 prune_segment(segmentsx,segmentx2);
837 }
838
839 npruned++;
840 }
841
842 endloop:
843
844 if(!((i+1)%10000))
845 printf_middle("Pruning Short Segments: Segments=%"Pindex_t" Short=%"Pindex_t" Pruned=%"Pindex_t,i+1,nshort,npruned);
846 }
847
848 /* Unmap from memory / close the files */
849
850 #if !SLIM
851 nodesx->data=UnmapFile(nodesx->data);
852 segmentsx->data=UnmapFile(segmentsx->data);
853 waysx->data=UnmapFile(waysx->data);
854 #else
855 nodesx->fd=SlimUnmapFile(nodesx->fd);
856 segmentsx->fd=SlimUnmapFile(segmentsx->fd);
857 waysx->fd=SlimUnmapFile(waysx->fd);
858 #endif
859
860 /* Print the final message */
861
862 printf_last("Pruned Short Segments: Segments=%"Pindex_t" Short=%"Pindex_t" Pruned=%"Pindex_t,segmentsx->number,nshort,npruned);
863 }
864
865
866 /*++++++++++++++++++++++++++++++++++++++
867 Prune out any nodes from straight highways where the introduced error is smaller than a specified maximum.
868
869 NodesX *nodesx The set of nodes to use.
870
871 SegmentsX *segmentsx The set of segments to use.
872
873 WaysX *waysx The set of ways to use.
874
875 distance_t maximum The maximum error to introduce.
876 ++++++++++++++++++++++++++++++++++++++*/
877
878 void PruneStraightHighwayNodes(NodesX *nodesx,SegmentsX *segmentsx,WaysX *waysx,distance_t maximum)
879 {
880 index_t i;
881 index_t npruned=0;
882 index_t nalloc;
883 BitMask *checked;
884 index_t *nodes,*segments;
885 double *lats,*lons;
886 double maximumf;
887
888 if(nodesx->number==0 || segmentsx->number==0 || waysx->number==0)
889 return;
890
891 /* Print the start message */
892
893 printf_first("Pruning Straight Highway Nodes: Nodes=0 Pruned=0");
894
895 /* Map into memory / open the files */
896
897 #if !SLIM
898 nodesx->data=MapFile(nodesx->filename_tmp);
899 segmentsx->data=MapFileWriteable(segmentsx->filename_tmp);
900 waysx->data=MapFile(waysx->filename_tmp);
901 #else
902 nodesx->fd=SlimMapFile(nodesx->filename_tmp);
903 segmentsx->fd=SlimMapFileWriteable(segmentsx->filename_tmp);
904 waysx->fd=SlimMapFile(waysx->filename_tmp);
905
906 InvalidateNodeXCache(nodesx->cache);
907 InvalidateSegmentXCache(segmentsx->cache);
908 InvalidateWayXCache(waysx->cache);
909 #endif
910
911 checked=AllocBitMask(nodesx->number);
912
913 log_malloc(checked,LengthBitMask(nodesx->number)*sizeof(BitMask));
914
915 nodes =(index_t*)malloc_logassert((nalloc=1024)*sizeof(index_t));
916 segments=(index_t*)malloc_logassert( nalloc *sizeof(index_t));
917
918 lats=(double*)malloc_logassert(nalloc*sizeof(double));
919 lons=(double*)malloc_logassert(nalloc*sizeof(double));
920
921 /* Loop through the nodes and find stretches of simple highway for possible modification */
922
923 maximumf=distance_to_km(maximum);
924
925 for(i=0;i<nodesx->number;i++)
926 {
927 int lowerbounded=0,upperbounded=0;
928 index_t lower=nalloc/2,current=nalloc/2,upper=nalloc/2;
929
930 if(IsBitSet(checked,i))
931 goto endloop;
932
933 if(segmentsx->firstnode[i]==NO_SEGMENT)
934 goto endloop;
935
936 /* Find all connected nodes */
937
938 nodes[current]=i;
939
940 do
941 {
942 index_t node1=NO_NODE,node2=NO_NODE;
943 index_t segment1=NO_SEGMENT,segment2=NO_SEGMENT;
944 index_t way1=NO_WAY,way2=NO_WAY;
945 int segcount=0;
946 NodeX *nodex;
947
948 /* Get the node data */
949
950 nodex=LookupNodeX(nodesx,nodes[current],1);
951
952 lats[current]=latlong_to_radians(nodex->latitude);
953 lons[current]=latlong_to_radians(nodex->longitude);
954
955 /* Count the segments at the node if not forced to be an end node */
956
957 if(IsBitSet(checked,nodes[current]))
958 ;
959 else if(nodex->flags&NODE_MINIRNDBT)
960 ;
961 else if(nodex->flags&NODE_TURNRSTRCT2 || nodex->flags&NODE_TURNRSTRCT)
962 ;
963 else
964 {
965 SegmentX *segmentx;
966
967 /* Count the segments connected to the node */
968
969 segmentx=FirstSegmentX(segmentsx,nodes[current],3);
970
971 while(segmentx)
972 {
973 segcount++;
974
975 if(node1==NO_NODE)
976 {
977 segment1=IndexSegmentX(segmentsx,segmentx);
978 node1=OtherNode(segmentx,nodes[current]);
979 way1=segmentx->way;
980 }
981 else if(node2==NO_NODE)
982 {
983 segment2=IndexSegmentX(segmentsx,segmentx);
984 node2=OtherNode(segmentx,nodes[current]);
985 way2=segmentx->way;
986 }
987 else
988 break;
989
990 segmentx=NextSegmentX(segmentsx,segmentx,nodes[current]);
991 }
992 }
993
994 /* Check if allowed due to one-way properties */
995
996 if(segcount==2)
997 {
998 SegmentX *segmentx1,*segmentx2;
999
1000 segmentx1=LookupSegmentX(segmentsx,segment1,1);
1001 segmentx2=LookupSegmentX(segmentsx,segment2,2);
1002
1003 if(!IsOneway(segmentx1) && !IsOneway(segmentx2))
1004 ;
1005 else if(IsOneway(segmentx1) && IsOneway(segmentx2))
1006 {
1007 if(IsOnewayTo(segmentx1,nodes[current]) && !IsOnewayFrom(segmentx2,nodes[current])) /* S1 is one-way but S2 doesn't continue */
1008 segcount=0;
1009
1010 if(IsOnewayFrom(segmentx1,nodes[current]) && !IsOnewayTo(segmentx2,nodes[current])) /* S1 is one-way but S2 doesn't continue */
1011 segcount=0;
1012 }
1013 else
1014 segcount=0;
1015 }
1016
1017 /* Check if allowed due to highway properties and node restrictions */
1018
1019 if(segcount==2)
1020 {
1021 WayX *wayx1,*wayx2;
1022
1023 wayx1=LookupWayX(waysx,way1,1);
1024 wayx2=LookupWayX(waysx,way2,2);
1025
1026 if(WaysCompare(&wayx1->way,&wayx2->way))
1027 segcount=0;
1028
1029 if(wayx1->way.name!=wayx2->way.name)
1030 segcount=0;
1031
1032 if((nodex->allow&wayx1->way.allow)!=wayx1->way.allow)
1033 segcount=0;
1034
1035 if((nodex->allow&wayx2->way.allow)!=wayx2->way.allow)
1036 segcount=0;
1037 }
1038
1039 /* Update the lists */
1040
1041 if(segcount==2)
1042 {
1043 /* Make space in the lists */
1044
1045 if(upper==(nalloc-1))
1046 {
1047 nodes =(index_t*)realloc_logassert(nodes ,(nalloc+=1024)*sizeof(index_t));
1048 segments=(index_t*)realloc_logassert(segments, nalloc *sizeof(index_t));
1049
1050 lats=(double*)realloc_logassert(lats,nalloc*sizeof(double));
1051 lons=(double*)realloc_logassert(lons,nalloc*sizeof(double));
1052 }
1053
1054 if(lower==0) /* move everything up by one */
1055 {
1056 memmove(nodes+1 ,nodes ,(1+upper-lower)*sizeof(index_t));
1057 memmove(segments+1,segments,(1+upper-lower)*sizeof(index_t));
1058
1059 memmove(lats+1,lats,(1+upper-lower)*sizeof(double));
1060 memmove(lons+1,lons,(1+upper-lower)*sizeof(double));
1061
1062 current++;
1063 lower++;
1064 upper++;
1065 }
1066
1067 if(lower==upper) /* first */
1068 {
1069 lower--;
1070
1071 nodes[lower]=node1;
1072 segments[lower]=segment1;
1073
1074 upper++;
1075
1076 nodes[upper]=node2;
1077 segments[upper-1]=segment2;
1078 segments[upper]=NO_SEGMENT;
1079
1080 current--;
1081 }
1082 else if(current==lower)
1083 {
1084 lower--;
1085
1086 if(nodes[current+1]==node2)
1087 {
1088 nodes[lower]=node1;
1089 segments[lower]=segment1;
1090 }
1091 else /* if(nodes[current+1]==node1) */
1092 {
1093 nodes[lower]=node2;
1094 segments[lower]=segment2;
1095 }
1096
1097 current--;
1098 }
1099 else /* if(current==upper) */
1100 {
1101 upper++;
1102
1103 if(nodes[current-1]==node2)
1104 {
1105 nodes[upper]=node1;
1106 segments[upper-1]=segment1;
1107 }
1108 else /* if(nodes[current-1]==node1) */
1109 {
1110 nodes[upper]=node2;
1111 segments[upper-1]=segment2;
1112 }
1113
1114 segments[upper]=NO_SEGMENT;
1115
1116 current++;
1117 }
1118
1119 if(nodes[upper]==nodes[lower])
1120 {
1121 if(!lowerbounded && !upperbounded)
1122 {
1123 nodex=LookupNodeX(nodesx,nodes[lower],1);
1124
1125 lats[lower]=latlong_to_radians(nodex->latitude);
1126 lons[lower]=latlong_to_radians(nodex->longitude);
1127 }
1128
1129 lats[upper]=lats[lower];
1130 lons[upper]=lons[lower];
1131
1132 lowerbounded=1;
1133 upperbounded=1;
1134 }
1135 }
1136 else /* if(segment!=2) */
1137 {
1138 if(current==upper)
1139 upperbounded=1;
1140
1141 if(current==lower)
1142 {
1143 lowerbounded=1;
1144 current=upper;
1145 }
1146 }
1147 }
1148 while(!(lowerbounded && upperbounded));
1149
1150 /* Mark the nodes */
1151
1152 for(current=lower;current<=upper;current++)
1153 SetBit(checked,nodes[current]);
1154
1155 /* Check for straight highway */
1156
1157 for(;lower<(upper-1);lower++)
1158 {
1159 for(current=upper;current>(lower+1);current--)
1160 {
1161 SegmentX *segmentx;
1162 distance_t dist=0;
1163 double dist1,dist2,dist3,distp;
1164 index_t c;
1165
1166 dist3=distance(lats[lower],lons[lower],lats[current],lons[current]);
1167
1168 for(c=lower+1;c<current;c++)
1169 {
1170 dist1=distance(lats[lower] ,lons[lower] ,lats[c],lons[c]);
1171 dist2=distance(lats[current],lons[current],lats[c],lons[c]);
1172
1173 /* Use law of cosines (assume flat Earth) */
1174
1175 if(dist3==0)
1176 distp=dist1; /* == dist2 */
1177 else if((dist1+dist2)<dist3)
1178 distp=0;
1179 else
1180 {
1181 double dist3a=(dist1*dist1-dist2*dist2+dist3*dist3)/(2.0*dist3);
1182 double dist3b=dist3-dist3a;
1183
1184 if(dist3a>=0 && dist3b>=0)
1185 distp=sqrt(dist1*dist1-dist3a*dist3a);
1186 else if(dist3a>0)
1187 distp=dist2;
1188 else /* if(dist3b>0) */
1189 distp=dist1;
1190 }
1191
1192 if(distp>maximumf) /* gone too far */
1193 break;
1194 }
1195
1196 if(c<current) /* not finished */
1197 continue;
1198
1199 /* Delete some segments and shift along */
1200
1201 for(c=lower+1;c<current;c++)
1202 {
1203 segmentx=LookupSegmentX(segmentsx,segments[c],1);
1204
1205 dist+=DISTANCE(segmentx->distance);
1206
1207 prune_segment(segmentsx,segmentx);
1208
1209 npruned++;
1210 }
1211
1212 segmentx=LookupSegmentX(segmentsx,segments[lower],1);
1213
1214 if(nodes[lower]==nodes[current]) /* loop; all within maximum distance */
1215 {
1216 prune_segment(segmentsx,segmentx);
1217
1218 npruned++;
1219 }
1220 else
1221 {
1222 segmentx->distance+=dist;
1223
1224 if(segmentx->node1==nodes[lower])
1225 modify_segment(segmentsx,segmentx,nodes[lower],nodes[current]);
1226 else /* if(segmentx->node2==nodes[lower]) */
1227 modify_segment(segmentsx,segmentx,nodes[current],nodes[lower]);
1228 }
1229
1230 lower=current-1;
1231 break;
1232 }
1233 }
1234
1235 endloop:
1236
1237 if(!((i+1)%10000))
1238 printf_middle("Pruning Straight Highway Nodes: Nodes=%"Pindex_t" Pruned=%"Pindex_t,i+1,npruned);
1239 }
1240
1241 /* Unmap from memory / close the files */
1242
1243 log_free(checked);
1244 free(checked);
1245
1246 free(nodes);
1247 free(segments);
1248
1249 free(lats);
1250 free(lons);
1251
1252 #if !SLIM
1253 nodesx->data=UnmapFile(nodesx->data);
1254 segmentsx->data=UnmapFile(segmentsx->data);
1255 waysx->data=UnmapFile(waysx->data);
1256 #else
1257 nodesx->fd=SlimUnmapFile(nodesx->fd);
1258 segmentsx->fd=SlimUnmapFile(segmentsx->fd);
1259 waysx->fd=SlimUnmapFile(waysx->fd);
1260 #endif
1261
1262 /* Print the final message */
1263
1264 printf_last("Pruned Straight Highway Nodes: Nodes=%"Pindex_t" Pruned=%"Pindex_t,nodesx->number,npruned);
1265 }
1266
1267
1268 /*++++++++++++++++++++++++++++++++++++++
1269 Prune a segment; unused nodes and ways will get marked for pruning later.
1270
1271 SegmentsX *segmentsx The set of segments to use.
1272
1273 SegmentX *segmentx The segment to be pruned.
1274 ++++++++++++++++++++++++++++++++++++++*/
1275
1276 static void prune_segment(SegmentsX *segmentsx,SegmentX *segmentx)
1277 {
1278 unlink_segment_node1_refs(segmentsx,segmentx);
1279
1280 unlink_segment_node2_refs(segmentsx,segmentx);
1281
1282 segmentx->node1=NO_NODE;
1283 segmentx->node2=NO_NODE;
1284 segmentx->next2=NO_SEGMENT;
1285
1286 PutBackSegmentX(segmentsx,segmentx);
1287 }
1288
1289
1290 /*++++++++++++++++++++++++++++++++++++++
1291 Modify a segment's nodes; unused nodes will get marked for pruning later.
1292
1293 SegmentsX *segmentsx The set of segments to use.
1294
1295 SegmentX *segmentx The segment to be modified.
1296
1297 index_t newnode1 The new value of node1.
1298
1299 index_t newnode2 The new value of node2.
1300 ++++++++++++++++++++++++++++++++++++++*/
1301
1302 static void modify_segment(SegmentsX *segmentsx,SegmentX *segmentx,index_t newnode1,index_t newnode2)
1303 {
1304 index_t thissegment=IndexSegmentX(segmentsx,segmentx);
1305
1306 if(newnode1>newnode2) /* rotate the segment around */
1307 {
1308 index_t temp;
1309
1310 if(segmentx->distance&(ONEWAY_2TO1|ONEWAY_1TO2))
1311 segmentx->distance^=ONEWAY_2TO1|ONEWAY_1TO2;
1312
1313 temp=newnode1;
1314 newnode1=newnode2;
1315 newnode2=temp;
1316 }
1317
1318 if(newnode1!=segmentx->node1)
1319 unlink_segment_node1_refs(segmentsx,segmentx);
1320
1321 if(newnode2!=segmentx->node2)
1322 unlink_segment_node2_refs(segmentsx,segmentx);
1323
1324 if(newnode1!=segmentx->node1) /* only modify it if the node has changed */
1325 {
1326 segmentx->node1=newnode1;
1327
1328 segmentsx->next1[thissegment]=segmentsx->firstnode[newnode1];
1329 segmentsx->firstnode[newnode1]=thissegment;
1330 }
1331
1332 if(newnode2!=segmentx->node2) /* only modify it if the node has changed */
1333 {
1334 segmentx->node2=newnode2;
1335
1336 segmentx->next2=segmentsx->firstnode[newnode2];
1337 segmentsx->firstnode[newnode2]=thissegment;
1338 }
1339
1340 PutBackSegmentX(segmentsx,segmentx);
1341 }
1342
1343
1344 /*++++++++++++++++++++++++++++++++++++++
1345 Unlink a node1 from a segment by modifying the linked list type arrangement of node references.
1346
1347 SegmentsX *segmentsx The set of segments to use.
1348
1349 SegmentX *segmentx The segment to be modified.
1350 ++++++++++++++++++++++++++++++++++++++*/
1351
1352 static void unlink_segment_node1_refs(SegmentsX *segmentsx,SegmentX *segmentx)
1353 {
1354 index_t segment,thissegment;
1355
1356 thissegment=IndexSegmentX(segmentsx,segmentx);
1357
1358 segment=segmentsx->firstnode[segmentx->node1];
1359
1360 if(segment==thissegment)
1361 segmentsx->firstnode[segmentx->node1]=segmentsx->next1[thissegment];
1362 else
1363 {
1364 do
1365 {
1366 index_t nextsegment;
1367 SegmentX *segx=LookupSegmentX(segmentsx,segment,4);
1368
1369 if(segx->node1==segmentx->node1)
1370 {
1371 nextsegment=segmentsx->next1[segment];
1372
1373 if(nextsegment==thissegment)
1374 segmentsx->next1[segment]=segmentsx->next1[thissegment];
1375 }
1376 else /* if(segx->node2==segmentx->node1) */
1377 {
1378 nextsegment=segx->next2;
1379
1380 if(nextsegment==thissegment)
1381 {
1382 segx->next2=segmentsx->next1[thissegment];
1383
1384 PutBackSegmentX(segmentsx,segx);
1385 }
1386 }
1387
1388 segment=nextsegment;
1389 }
1390 while(segment!=thissegment && segment!=NO_SEGMENT);
1391 }
1392 }
1393
1394
1395 /*++++++++++++++++++++++++++++++++++++++
1396 Unlink a node2 from a segment by modifying the linked list type arrangement of node references.
1397
1398 SegmentsX *segmentsx The set of segments to use.
1399
1400 SegmentX *segmentx The segment to be modified.
1401 ++++++++++++++++++++++++++++++++++++++*/
1402
1403 static void unlink_segment_node2_refs(SegmentsX *segmentsx,SegmentX *segmentx)
1404 {
1405 index_t segment,thissegment;
1406
1407 thissegment=IndexSegmentX(segmentsx,segmentx);
1408
1409 segment=segmentsx->firstnode[segmentx->node2];
1410
1411 if(segment==thissegment)
1412 segmentsx->firstnode[segmentx->node2]=segmentx->next2;
1413 else
1414 {
1415 do
1416 {
1417 index_t nextsegment;
1418 SegmentX *segx=LookupSegmentX(segmentsx,segment,4);
1419
1420 if(segx->node1==segmentx->node2)
1421 {
1422 nextsegment=segmentsx->next1[segment];
1423
1424 if(nextsegment==thissegment)
1425 segmentsx->next1[segment]=segmentx->next2;
1426 }
1427 else /* if(segx->node2==segmentx->node2) */
1428 {
1429 nextsegment=segx->next2;
1430
1431 if(nextsegment==thissegment)
1432 {
1433 segx->next2=segmentx->next2;
1434
1435 PutBackSegmentX(segmentsx,segx);
1436 }
1437 }
1438
1439 segment=nextsegment;
1440 }
1441 while(segment!=thissegment && segment!=NO_SEGMENT);
1442 }
1443 }
1444
1445
1446 /*++++++++++++++++++++++++++++++++++++++
1447 Calculate the distance between two locations.
1448
1449 double distance Returns the distance between the locations.
1450
1451 double lat1 The latitude of the first location.
1452
1453 double lon1 The longitude of the first location.
1454
1455 double lat2 The latitude of the second location.
1456
1457 double lon2 The longitude of the second location.
1458 ++++++++++++++++++++++++++++++++++++++*/
1459
1460 static double distance(double lat1,double lon1,double lat2,double lon2)
1461 {
1462 double dlon = lon1 - lon2;
1463 double dlat = lat1 - lat2;
1464
1465 double a1,a2,a,sa,c,d;
1466
1467 if(dlon==0 && dlat==0)
1468 return 0;
1469
1470 a1 = sin (dlat / 2);
1471 a2 = sin (dlon / 2);
1472 a = (a1 * a1) + cos (lat1) * cos (lat2) * a2 * a2;
1473 sa = sqrt (a);
1474 if (sa <= 1.0)
1475 {c = 2 * asin (sa);}
1476 else
1477 {c = 2 * asin (1.0);}
1478 d = 6378.137 * c;
1479
1480 return(d);
1481 }