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 1417 - (show annotations) (download) (as text)
Sat Jun 22 18:55:34 2013 UTC (11 years, 9 months ago) by amb
File MIME type: text/x-csrc
File size: 41922 byte(s)
Use SlimMapFile(), SlimUnmapFile(), SlimFetch() and SlimReplace() [see earlier
log message] and also refactor the code so that the additional ways are placed
in a separate file using the Slim*() functions rather than being appended to the
open (or mapped) way file.

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