1 /*
2 * geoLines.c --
3 *
4 * This file defines functions that manage geolines and geolinearrays.
5 * See the user documentation for more information.
6 *
7 * Copyright (c) 2004 Gordon D. Carrie. All rights reserved.
8 *
9 * Licensed under the Open Software License version 2.1
10 *
11 * Please address questions and feedback to user0@tkgeomap.org
12 *
13 * @(#) $Id: geoLines.c,v 1.7 2005/08/29 16:01:52 tkgeomap Exp $
14 *
15 ********************************************
16 *
17 */
18
19 #include <stdlib.h>
20 #include <string.h>
21 #include <limits.h>
22 #include <math.h>
23
24 #include <tcl.h>
25 #include "geoLines.h"
26
27 /*
28 * Local function declaration.
29 */
30
31 static void lnInit _ANSI_ARGS_((struct GeoLn *lnPtr));
32
33
34 /*
35 *----------------------------------------------------------------------
36 *
37 * GeoLnCreate --
38 *
39 * This procedure creates and initializes a new geoline.
40 *
41 * Results:
42 * See the user documentation.
43 *
44 * Side effects:
45 * See the user documentation.
46 *
47 *----------------------------------------------------------------------
48 */
49
50 GeoLn
GeoLnCreate(nPtsMax)51 GeoLnCreate(nPtsMax)
52 unsigned nPtsMax;
53 {
54 struct GeoLn *lnPtr;
55
56 lnPtr = (GeoLn)CKALLOC(sizeof(*lnPtr));
57 lnInit(lnPtr);
58 if (nPtsMax == 0) {
59 return lnPtr;
60 }
61 lnPtr->pts = (GeoPt *)CKALLOC(nPtsMax * sizeof(GeoPt));
62 lnPtr->nPtsMax = nPtsMax;
63 return lnPtr;
64 }
65
66 /*
67 *----------------------------------------------------------------------
68 *
69 * lnInit --
70 *
71 * This procedure initializes a new geoline.
72 *
73 * Results:
74 * None.
75 *
76 * Side effects:
77 * The contents of a GeoLn structure are set. The previous contents,
78 * which are assumed to be garbage, are ignored.
79 *
80 *----------------------------------------------------------------------
81 */
82
83 static void
lnInit(lnPtr)84 lnInit(lnPtr)
85 struct GeoLn *lnPtr;
86 {
87 /* Initialize a new Line struct. */
88
89 lnPtr->nPts = lnPtr->nPtsMax = 0;
90 lnPtr->lonMax = lnPtr->latMax = -INT_MAX;
91 lnPtr->lonMin = lnPtr->latMin = INT_MAX;
92 lnPtr->pts = NULL;
93 }
94
95 /*
96 *----------------------------------------------------------------------
97 *
98 * GeoLnClear --
99 *
100 * This procedure set the number of points in a geoline to zero without
101 * releasing its internal storage.
102 *
103 * Results:
104 * None.
105 *
106 * Side effects:
107 * See the user documentation.
108 *
109 *----------------------------------------------------------------------
110 */
111
112 void
GeoLnClear(lnPtr)113 GeoLnClear(lnPtr)
114 struct GeoLn *lnPtr;
115 {
116 lnPtr->nPts = 0;
117 lnPtr->lonMax = lnPtr->latMax = -INT_MAX;
118 lnPtr->lonMin = lnPtr->latMin = INT_MAX;
119 }
120
121 /*
122 *----------------------------------------------------------------------
123 *
124 * GeoLnDestroy --
125 *
126 * This procedure frees internal storage in a GeoLn structure, and
127 * frees the structure.
128 *
129 * Results:
130 * None.
131 *
132 * Side effects:
133 * See the user documentation.
134 *
135 *----------------------------------------------------------------------
136 */
137
138 void
GeoLnDestroy(lnPtr)139 GeoLnDestroy(lnPtr)
140 struct GeoLn *lnPtr;
141 {
142 if (!lnPtr) {
143 return;
144 }
145 CKFREE((char *)lnPtr->pts);
146 CKFREE((char *)lnPtr);
147 }
148
149 /*
150 *----------------------------------------------------------------------
151 *
152 * GeoLnSetAlloc --
153 *
154 * This procedure changes the number of points a geoline can store.
155 *
156 * Results:
157 * None.
158 *
159 * Side effects:
160 * This procedure reallocates a GeoLn structure's internal storage.
161 * If the number of points is reduced, the geoline's maxima and minima
162 * are recomputed.
163 *
164 *----------------------------------------------------------------------
165 */
166
167 void
GeoLnSetAlloc(lnPtr,nPtsMax)168 GeoLnSetAlloc(lnPtr, nPtsMax)
169 struct GeoLn *lnPtr;
170 unsigned nPtsMax;
171 {
172 if (lnPtr->nPtsMax == nPtsMax) {
173 return;
174 }
175 if (nPtsMax == 0) {
176 CKFREE((char *)lnPtr->pts);
177 lnInit(lnPtr);
178 return;
179 }
180 lnPtr->pts = (GeoPt *)CKREALLOC((char *)lnPtr->pts,
181 nPtsMax * sizeof(GeoPt));
182 lnPtr->nPtsMax = nPtsMax;
183 if (lnPtr->nPts > lnPtr->nPtsMax) {
184 /*
185 * If shrinking line, reset nPts and recompute maxima and minima.
186 */
187
188 GeoPt *p, *pe;
189 lnPtr->nPts = lnPtr->nPtsMax;
190 for (p = lnPtr->pts, pe = lnPtr->pts + lnPtr->nPts; p < pe; p++) {
191 if (GeoPtIsSomewhere(*p)) {
192 lnPtr->latMax
193 = (lnPtr->latMax > p->lat) ? lnPtr->latMax : p->lat;
194 lnPtr->lonMax
195 = (lnPtr->lonMax > p->lon) ? lnPtr->lonMax : p->lon;
196 lnPtr->latMin
197 = (lnPtr->latMin < p->lat) ? lnPtr->latMin : p->lat;
198 lnPtr->lonMin
199 = (lnPtr->lonMin < p->lon) ? lnPtr->lonMin : p->lon;
200 }
201 }
202 }
203 }
204
205 /*
206 *----------------------------------------------------------------------
207 *
208 * GeoLnAddPt --
209 *
210 * This procedure adds a point to a geoline.
211 *
212 * Results:
213 * None.
214 *
215 * Side effects:
216 * The geoline's allocation might increase. It's maxima and minima will
217 * be adjusted.
218 *
219 *----------------------------------------------------------------------
220 */
221
222 void
GeoLnAddPt(p,lnPtr)223 GeoLnAddPt(p, lnPtr)
224 GeoPt p;
225 struct GeoLn *lnPtr;
226 {
227 if (lnPtr->nPts + 1 > lnPtr->nPtsMax) {
228 GeoLnSetAlloc(lnPtr, ((lnPtr->nPtsMax + 4) * 5) / 4);
229 }
230 if (GeoPtIsSomewhere(p)) {
231 lnPtr->latMax = (lnPtr->latMax > p.lat) ? lnPtr->latMax : p.lat;
232 lnPtr->lonMax = (lnPtr->lonMax > p.lon) ? lnPtr->lonMax : p.lon;
233 lnPtr->latMin = (lnPtr->latMin < p.lat) ? lnPtr->latMin : p.lat;
234 lnPtr->lonMin = (lnPtr->lonMin < p.lon) ? lnPtr->lonMin : p.lon;
235 }
236 lnPtr->pts[lnPtr->nPts] = p;
237 lnPtr->nPts++;
238 return;
239 }
240
241 /*
242 *----------------------------------------------------------------------
243 *
244 * GeoLnGetPt --
245 *
246 * This procedure retrieve a geopoint from a geoline.
247 *
248 * Results:
249 * See the user documentation.
250 *
251 * Side effects:
252 * None.
253 *
254 *----------------------------------------------------------------------
255 */
256
257 GeoPt
GeoLnGetPt(lnPtr,n)258 GeoLnGetPt(lnPtr, n)
259 struct GeoLn *lnPtr;
260 unsigned n;
261 {
262 /* Return MultiPt n from lnPtr, or a bogus location if out of bounds. */
263 return (n < lnPtr->nPts) ? *(lnPtr->pts + n) : GeoPtNowhere();
264 }
265
266 /*
267 *----------------------------------------------------------------------
268 *
269 * GeoLnCtr --
270 *
271 * This procedure computes the average position of the points in a
272 * geoline in 3-dimensional space.
273 *
274 * Results:
275 * See the user documentation.
276 *
277 * Side effects:
278 * None.
279 *
280 *----------------------------------------------------------------------
281 */
282
283 CartPt
GeoLnCtr(lnPtr)284 GeoLnCtr(lnPtr)
285 struct GeoLn *lnPtr;
286 {
287 CartPt ctr = {0.0, 0.0, 0.0}; /* Return value */
288 CartPt cpt;
289 GeoPt *p, *end; /* Loop index */
290
291 for (p = lnPtr->pts, end = p + lnPtr->nPts; p < end; p++) {
292 cpt = LatLonToCart(*p);
293 ctr.x += cpt.x;
294 ctr.y += cpt.y;
295 ctr.z += cpt.z;
296 }
297 ctr.x /= lnPtr->nPts;
298 ctr.y /= lnPtr->nPts;
299 ctr.z /= lnPtr->nPts;
300 return ctr;
301 }
302
303 /*
304 *----------------------------------------------------------------------
305 *
306 * GeoLnContainGeoPt --
307 *
308 * This procedure determines if a geoline enclosed a geopoint.
309 *
310 * Results:
311 * See the user documentation.
312 *
313 * Side effects:
314 * None.
315 *
316 *----------------------------------------------------------------------
317 */
318
319 int
GeoLnContainGeoPt(geoPt,lnPtr)320 GeoLnContainGeoPt(geoPt, lnPtr)
321 GeoPt geoPt;
322 struct GeoLn *lnPtr;
323 {
324 int mrdx; /* Number of times line crosses meridian containing geoPt */
325 int lnx; /* Number of times line crosses line from geoPt to North pole */
326 GeoPt *p0, *p1, *end;
327
328 /*
329 * Loop through segments in lnPtr, counting number of times lnPtr
330 * crosses meridian, and number of crossings between geoPt and North pole
331 */
332
333 for (mrdx = 0, lnx = 0,
334 p0 = lnPtr->pts + lnPtr->nPts - 1,
335 p1 = lnPtr->pts,
336 end = lnPtr->pts + lnPtr->nPts;
337 p1 < end;
338 p0 = p1++) {
339 /*
340 * Determine if segment defined by p0--p1 straddles meridian
341 * containing geoPt, or is on boundary. Do not count segments on
342 * boundary as more than one crossing
343 */
344
345 if (LonBtwn1(geoPt.lon, p1->lon, p0->lon)) {
346 double lat0 = AngleToDeg(p0->lat);
347 double lon0 = AngleToDeg(p0->lon);
348 double lat1 = AngleToDeg(p1->lat);
349 double lon1 = AngleToDeg(p1->lon);
350 double xlat; /* Latitude of segment crossing */
351
352 mrdx++;
353 xlat = lat0 + (AngleToDeg(geoPt.lon) - lon0)
354 * (lat1 - lat0) / (lon1 - lon0);
355 if (LatCmp(AngleFmDeg(xlat), geoPt.lat) == North) {
356 lnx = !lnx;
357 }
358
359 }
360 }
361
362 if (mrdx % 2 == 1) {
363 /*
364 * Odd number of meridian crossings => region contains a pole */
365 /* Assume pole is that of hemisphere containing lnPtr's mean
366 */
367
368 CartPt ctr = GeoLnCtr(lnPtr);
369 if (ctr.z > 0.0) {
370 /*
371 * Region contains North Pole
372 */
373
374 return !lnx;
375 }
376
377 }
378 return lnx;
379 }
380
381 /*
382 *----------------------------------------------------------------------
383 *
384 * GeoLnArrCreate --
385 *
386 * This procedure creates and initializes an new geolinearray.
387 *
388 * Results:
389 * See the user documentation.
390 *
391 * Side effects:
392 * See the user documentation.
393 *
394 *----------------------------------------------------------------------
395 */
396
397 GeoLnArr
GeoLnArrCreate(nLinesMax)398 GeoLnArrCreate(nLinesMax)
399 unsigned nLinesMax;
400 {
401 struct GeoLnArr *lnArrPtr;
402 int n;
403
404 lnArrPtr = (GeoLnArr)CKALLOC(sizeof(*lnArrPtr));
405 lnArrPtr->descr = NULL;
406 lnArrPtr->lines = NULL;
407 GeoLnArrSetDescr(lnArrPtr, "");
408 lnArrPtr->nLines = lnArrPtr->nLinesMax = 0;
409 lnArrPtr->nPts = lnArrPtr->nMax = 0;
410 lnArrPtr->lonMax = lnArrPtr->latMax = -INT_MAX;
411 lnArrPtr->lonMin = lnArrPtr->latMin = INT_MAX;
412 lnArrPtr->lines = NULL;
413 if (nLinesMax == 0) {
414 return lnArrPtr;
415 }
416 lnArrPtr->lines = (struct GeoLn **)CKALLOC(nLinesMax
417 * sizeof(struct GeoLn *));
418 lnArrPtr->nLinesMax = nLinesMax;
419 for (n = 0; n < nLinesMax; n++) {
420 lnArrPtr->lines[n] = NULL;
421 }
422 return lnArrPtr;
423 }
424
425 /*
426 *----------------------------------------------------------------------
427 *
428 * GeoLnArrSetDescr --
429 *
430 * This procedure sets the descriptor for a geolinearray.
431 *
432 * Results:
433 * None.
434 *
435 * Side effects:
436 * See the user documentation.
437 *
438 *----------------------------------------------------------------------
439 */
440
441 void
GeoLnArrSetDescr(lnArrPtr,descr)442 GeoLnArrSetDescr(lnArrPtr, descr)
443 struct GeoLnArr *lnArrPtr;
444 CONST char *descr;
445 {
446 lnArrPtr->descr = CKREALLOC(lnArrPtr->descr, strlen(descr) + 1);
447 strcpy(lnArrPtr->descr, descr);
448 }
449
450 /*
451 *----------------------------------------------------------------------
452 *
453 * GeoLnArrSetAlloc --
454 *
455 * This procedure changes the number of geolines a geolinearray can store.
456 *
457 * Results:
458 * None.
459 *
460 * Side effects:
461 * See the user documentation.
462 *
463 *----------------------------------------------------------------------
464 */
465
466 void
GeoLnArrSetAlloc(lnArrPtr,nLinesMax)467 GeoLnArrSetAlloc(lnArrPtr, nLinesMax)
468 struct GeoLnArr *lnArrPtr;
469 unsigned nLinesMax;
470 {
471 int n; /* Loop index */
472
473 if (lnArrPtr->nLinesMax == nLinesMax) {
474 return;
475 }
476 for (n = nLinesMax; n < lnArrPtr->nLinesMax; n++) {
477 /*
478 * Free excess lines
479 */
480
481 GeoLnDestroy(lnArrPtr->lines[n]);
482 }
483 lnArrPtr->lines = (GeoLn *)CKREALLOC((char *)lnArrPtr->lines,
484 nLinesMax * sizeof(GeoLn));
485 lnArrPtr->nLinesMax = nLinesMax;
486 for (n = lnArrPtr->nLines; n < lnArrPtr->nLinesMax; n++) {
487 /*
488 * Initialize additional lines.
489 */
490
491 lnArrPtr->lines[n] = NULL;
492 }
493
494 }
495
496 /*
497 *----------------------------------------------------------------------
498 *
499 * GeoLnArrAddLine --
500 *
501 * This procedure copies a geoline onto the end of a geolinearray.
502 *
503 * Results:
504 * See the user documentation.
505 *
506 * Side effects:
507 * See the user documentation.
508 *
509 *----------------------------------------------------------------------
510 */
511
512 int
GeoLnArrAddLine(ln,lnArrPtr)513 GeoLnArrAddLine(ln, lnArrPtr)
514 GeoLn ln;
515 struct GeoLnArr *lnArrPtr;
516 {
517 int nLines; /* Initial number of lines in lnArrPtr */
518
519 nLines = lnArrPtr->nLines;
520 if (nLines + 1 > lnArrPtr->nLinesMax) {
521 /*
522 * If lnArrPtr needs more space in its allocation, grow it by 25%
523 */
524
525 GeoLnArrSetAlloc(lnArrPtr, ((lnArrPtr->nLinesMax + 4) * 5) / 4);
526 }
527 if ( !(lnArrPtr->lines[nLines] = GeoLnCreate(ln->nPts)) ) {
528 return 0;
529 }
530 lnArrPtr->nPts += ln->nPts;
531 lnArrPtr->nMax = (lnArrPtr->nMax > ln->nPts) ? lnArrPtr->nMax : ln->nPts;
532 lnArrPtr->latMax
533 = (lnArrPtr->latMax > ln->latMax) ? lnArrPtr->latMax : ln->latMax;
534 lnArrPtr->lonMax
535 = (lnArrPtr->lonMax > ln->lonMax) ? lnArrPtr->lonMax : ln->lonMax;
536 lnArrPtr->latMin
537 = (lnArrPtr->latMin < ln->latMin) ? lnArrPtr->latMin : ln->latMin;
538 lnArrPtr->lonMin
539 = (lnArrPtr->lonMin < ln->lonMin) ? lnArrPtr->lonMin : ln->lonMin;
540 memcpy(lnArrPtr->lines[nLines]->pts, ln->pts, ln->nPts * sizeof(GeoPt));
541 lnArrPtr->lines[nLines]->nPts = ln->nPts;
542 lnArrPtr->lines[nLines]->lonMax = ln->lonMax;
543 lnArrPtr->lines[nLines]->lonMin = ln->lonMin;
544 lnArrPtr->lines[nLines]->latMax = ln->latMax;
545 lnArrPtr->lines[nLines]->latMin = ln->latMin;
546 lnArrPtr->nLines++;
547 return 1;
548 }
549
550 /*
551 *----------------------------------------------------------------------
552 *
553 * GeoLnArrPutLine --
554 *
555 * This procedure gives a geoline to a geolinearray.
556 *
557 * Results:
558 * See the user documentation.
559 *
560 * Side effects:
561 * See the user documentation.
562 *
563 *----------------------------------------------------------------------
564 */
565
566 int
GeoLnArrPutLine(ln,lnArrPtr)567 GeoLnArrPutLine(ln, lnArrPtr)
568 GeoLn ln;
569 struct GeoLnArr *lnArrPtr;
570 {
571 int nLines; /* Initial number of lines in lnArrPtr */
572
573 nLines = lnArrPtr->nLines;
574 if (nLines + 1 > lnArrPtr->nLinesMax) {
575 /*
576 * If lnArrPtr needs more space in its allocation, grow it by 25%
577 */
578
579 GeoLnArrSetAlloc(lnArrPtr, ((lnArrPtr->nLinesMax + 4) * 5) / 4);
580 }
581 lnArrPtr->nPts += ln->nPts;
582 lnArrPtr->nMax = (lnArrPtr->nMax > ln->nPts) ? lnArrPtr->nMax : ln->nPts;
583 lnArrPtr->latMax
584 = (lnArrPtr->latMax > ln->latMax) ? lnArrPtr->latMax : ln->latMax;
585 lnArrPtr->lonMax
586 = (lnArrPtr->lonMax > ln->lonMax) ? lnArrPtr->lonMax : ln->lonMax;
587 lnArrPtr->latMin
588 = (lnArrPtr->latMin < ln->latMin) ? lnArrPtr->latMin : ln->latMin;
589 lnArrPtr->lonMin
590 = (lnArrPtr->lonMin < ln->lonMin) ? lnArrPtr->lonMin : ln->lonMin;
591 lnArrPtr->lines[nLines] = ln;
592 lnArrPtr->nLines++;
593 return 1;
594 }
595
596 /*
597 *----------------------------------------------------------------------
598 *
599 * GeoLnArrGetDescr --
600 *
601 * This is a member access function.
602 *
603 * Results:
604 * See the user documentation.
605 *
606 * Side effects:
607 * None.
608 *
609 *----------------------------------------------------------------------
610 */
611
612 char *
GeoLnArrGetDescr(lnArrPtr)613 GeoLnArrGetDescr(lnArrPtr)
614 struct GeoLnArr *lnArrPtr;
615 {
616 return lnArrPtr->descr;
617 }
618
619 /*
620 *----------------------------------------------------------------------
621 *
622 * GeoLnArrGetLine --
623 *
624 * This procedure retrieves a geoline from a geolinearray.
625 *
626 * Results:
627 * See the user documentation.
628 *
629 * Side effects:
630 * None.
631 *
632 *----------------------------------------------------------------------
633 */
634
635 GeoLn
GeoLnArrGetLine(lnArrPtr,n)636 GeoLnArrGetLine(lnArrPtr, n)
637 struct GeoLnArr *lnArrPtr;
638 unsigned n;
639 {
640 return (n < lnArrPtr->nLines) ? lnArrPtr->lines[n] : NULL;
641 }
642
643 /*
644 *----------------------------------------------------------------------
645 *
646 * GeoLnArrFree --
647 *
648 * This procedure frees internal storage in a geolinearray.
649 *
650 * Results:
651 * None.
652 *
653 * Side effects:
654 * See the user documentation.
655 *
656 *----------------------------------------------------------------------
657 */
658
659 void
GeoLnArrFree(lnArrPtr)660 GeoLnArrFree(lnArrPtr)
661 struct GeoLnArr *lnArrPtr;
662 {
663 int n;
664
665 if ( !lnArrPtr ) {
666 return;
667 }
668 for (n = 0; n < lnArrPtr->nLines; n++) {
669 GeoLnDestroy(lnArrPtr->lines[n]);
670 }
671 CKFREE((char *)lnArrPtr->lines);
672 CKFREE((char *)lnArrPtr->descr);
673 }
674
675 /*
676 *----------------------------------------------------------------------
677 *
678 * GeoLnArrDestroy --
679 *
680 * This procedure frees internal storage in a geolinearray and frees the
681 * geolinearray.
682 *
683 * Results:
684 * None.
685 *
686 * Side effects:
687 * See the user documentation.
688 *
689 *----------------------------------------------------------------------
690 */
691
692 void
GeoLnArrDestroy(lnArrPtr)693 GeoLnArrDestroy(lnArrPtr)
694 struct GeoLnArr *lnArrPtr;
695 {
696 GeoLnArrFree(lnArrPtr);
697 CKFREE((char *)lnArrPtr);
698 }
699
700 /*
701 *----------------------------------------------------------------------
702 *
703 * GeoLnArrContainGeoPt --
704 *
705 * This procedure determines whether any of the geolines in a geolinearray
706 * contain a geopoint.
707 *
708 * Results:
709 * See the user documentation.
710 *
711 * Side effects:
712 * None.
713 *
714 *----------------------------------------------------------------------
715 */
716
717 int
GeoLnArrContainGeoPt(geoPt,lnArrPtr)718 GeoLnArrContainGeoPt(geoPt, lnArrPtr)
719 GeoPt geoPt;
720 struct GeoLnArr *lnArrPtr;
721 {
722 int n;
723
724 for (n = 0; n < lnArrPtr->nLines; n++) {
725 if ( GeoLnContainGeoPt(geoPt, lnArrPtr->lines[n]) ) {
726 return 1;
727 }
728 }
729 return 0;
730 }
731