/*
 * geoLines.c --
 *
 * 	This file defines functions that manage geolines and geolinearrays.
 * 	See the user documentation for more information.
 * 
 * Copyright (c) 2004 Gordon D. Carrie.  All rights reserved.
 * 
 * Licensed under the Open Software License version 2.1
 * 
 * Please address questions and feedback to user0@tkgeomap.org
 *
 * @(#) $Id: geoLines.c,v 1.7 2005/08/29 16:01:52 tkgeomap Exp $
 *
 ********************************************
 *
 */

#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <math.h>

#include <tcl.h>
#include "geoLines.h"

/*
 * Local function declaration.
 */

static void lnInit _ANSI_ARGS_((struct GeoLn *lnPtr));


/*
 *----------------------------------------------------------------------
 *
 * GeoLnCreate --
 *
 *	This procedure creates and initializes a new geoline.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

GeoLn
GeoLnCreate(nPtsMax)
    unsigned nPtsMax;
{
    struct GeoLn *lnPtr;

    lnPtr = (GeoLn)CKALLOC(sizeof(*lnPtr));
    lnInit(lnPtr);
    if (nPtsMax == 0) {
	return lnPtr;
    }
    lnPtr->pts = (GeoPt *)CKALLOC(nPtsMax * sizeof(GeoPt));
    lnPtr->nPtsMax = nPtsMax;
    return lnPtr;
}

/*
 *----------------------------------------------------------------------
 *
 * lnInit --
 *
 *	This procedure initializes a new geoline.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	The contents of a GeoLn structure are set.  The previous contents,
 *	which are assumed to be garbage, are ignored.
 *
 *----------------------------------------------------------------------
 */

static void
lnInit(lnPtr)
    struct GeoLn *lnPtr;
{
    /* Initialize a new Line struct. */

    lnPtr->nPts = lnPtr->nPtsMax = 0;
    lnPtr->lonMax = lnPtr->latMax = -INT_MAX;
    lnPtr->lonMin = lnPtr->latMin = INT_MAX;
    lnPtr->pts = NULL;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnClear --
 *
 *	This procedure set the number of points in a geoline to zero without
 *	releasing its internal storage.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnClear(lnPtr)
    struct GeoLn *lnPtr;
{
    lnPtr->nPts = 0;
    lnPtr->lonMax = lnPtr->latMax = -INT_MAX;
    lnPtr->lonMin = lnPtr->latMin = INT_MAX;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnDestroy --
 *
 *	This procedure frees internal storage in a GeoLn structure, and
 *	frees the structure.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnDestroy(lnPtr)
    struct GeoLn *lnPtr;
{
    if (!lnPtr) {
	return;
    }
    CKFREE((char *)lnPtr->pts);
    CKFREE((char *)lnPtr);
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnSetAlloc --
 *
 *	This procedure changes the number of points a geoline can store.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	This procedure reallocates a GeoLn structure's internal storage.
 *	If the number of points is reduced, the geoline's maxima and minima
 *	are recomputed.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnSetAlloc(lnPtr, nPtsMax)
    struct GeoLn *lnPtr;
    unsigned nPtsMax;
{
    if (lnPtr->nPtsMax == nPtsMax) {
	return;
    }
    if (nPtsMax == 0) {
	CKFREE((char *)lnPtr->pts);
	lnInit(lnPtr);
	return;
    }
    lnPtr->pts = (GeoPt *)CKREALLOC((char *)lnPtr->pts,
	    nPtsMax * sizeof(GeoPt));
    lnPtr->nPtsMax = nPtsMax;
    if (lnPtr->nPts > lnPtr->nPtsMax) {
	/*
	 * If shrinking line, reset nPts and recompute maxima and minima.
	 */

	GeoPt *p, *pe;
	lnPtr->nPts = lnPtr->nPtsMax;
	for (p = lnPtr->pts, pe = lnPtr->pts + lnPtr->nPts; p < pe; p++) {
	    if (GeoPtIsSomewhere(*p)) {
		lnPtr->latMax
		    = (lnPtr->latMax > p->lat) ? lnPtr->latMax : p->lat;
		lnPtr->lonMax
		    = (lnPtr->lonMax > p->lon) ? lnPtr->lonMax : p->lon;
		lnPtr->latMin
		    = (lnPtr->latMin < p->lat) ? lnPtr->latMin : p->lat;
		lnPtr->lonMin
		    = (lnPtr->lonMin < p->lon) ? lnPtr->lonMin : p->lon;
	    }
	}
    }
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnAddPt --
 *
 *	This procedure adds a point to a geoline.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	The geoline's allocation might increase.  It's maxima and minima will
 *	be adjusted.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnAddPt(p, lnPtr)
    GeoPt p;
    struct GeoLn *lnPtr;
{
    if (lnPtr->nPts + 1 > lnPtr->nPtsMax) {
	GeoLnSetAlloc(lnPtr, ((lnPtr->nPtsMax + 4) * 5) / 4);
    }
    if (GeoPtIsSomewhere(p)) {
	lnPtr->latMax = (lnPtr->latMax > p.lat) ? lnPtr->latMax : p.lat;
	lnPtr->lonMax = (lnPtr->lonMax > p.lon) ? lnPtr->lonMax : p.lon;
	lnPtr->latMin = (lnPtr->latMin < p.lat) ? lnPtr->latMin : p.lat;
	lnPtr->lonMin = (lnPtr->lonMin < p.lon) ? lnPtr->lonMin : p.lon;
    }
    lnPtr->pts[lnPtr->nPts] = p;
    lnPtr->nPts++;
    return;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnGetPt --
 *
 *	This procedure retrieve a geopoint from a geoline.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

GeoPt
GeoLnGetPt(lnPtr, n)
    struct GeoLn *lnPtr;
    unsigned n;
{
    /* Return MultiPt n from lnPtr, or a bogus location if out of bounds. */
    return (n < lnPtr->nPts) ? *(lnPtr->pts + n) : GeoPtNowhere();
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnCtr --
 *
 *	This procedure computes the average position of the points in a
 *	geoline in 3-dimensional space.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

CartPt
GeoLnCtr(lnPtr)
    struct GeoLn *lnPtr;
{
    CartPt ctr = {0.0, 0.0, 0.0};	/* Return value */
    CartPt cpt;
    GeoPt *p, *end;			/* Loop index */

    for (p = lnPtr->pts, end = p + lnPtr->nPts; p < end; p++) {
	cpt = LatLonToCart(*p);
	ctr.x += cpt.x;
	ctr.y += cpt.y;
	ctr.z += cpt.z;
    }
    ctr.x /= lnPtr->nPts;
    ctr.y /= lnPtr->nPts;
    ctr.z /= lnPtr->nPts;
    return ctr;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnContainGeoPt --
 *
 *	This procedure determines if a geoline enclosed a geopoint.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

int
GeoLnContainGeoPt(geoPt, lnPtr)
    GeoPt geoPt;
    struct GeoLn *lnPtr;
{
    int mrdx;	/* Number of times line crosses meridian containing geoPt */
    int lnx;	/* Number of times line crosses line from geoPt to North pole */
    GeoPt *p0, *p1, *end;

    /*
     * Loop through segments in lnPtr, counting number of times lnPtr
     * crosses meridian, and number of crossings between geoPt and North pole
     */

    for (mrdx = 0, lnx = 0,
	    p0 = lnPtr->pts + lnPtr->nPts - 1,
	    p1 = lnPtr->pts,
	    end = lnPtr->pts + lnPtr->nPts;
	    p1 < end;
	    p0 = p1++) {
	/*
	 * Determine if segment defined by p0--p1 straddles meridian
	 * containing geoPt, or is on boundary.  Do not count segments on
	 * boundary as more than one crossing
	 */

	if (LonBtwn1(geoPt.lon, p1->lon, p0->lon)) {
	    double lat0 = AngleToDeg(p0->lat);
	    double lon0 = AngleToDeg(p0->lon);
	    double lat1 = AngleToDeg(p1->lat);
	    double lon1 = AngleToDeg(p1->lon);
	    double xlat;		/* Latitude of segment crossing */

	    mrdx++;
	    xlat = lat0 + (AngleToDeg(geoPt.lon) - lon0)
		* (lat1 - lat0) / (lon1 - lon0);
	    if (LatCmp(AngleFmDeg(xlat), geoPt.lat) == North) {
		lnx = !lnx;
	    }

	}
    }

    if (mrdx % 2 == 1) {
	/*
	 * Odd number of meridian crossings => region contains a pole */
	/* Assume pole is that of hemisphere containing lnPtr's mean
	 */

	CartPt ctr = GeoLnCtr(lnPtr);
	if (ctr.z > 0.0) {
	    /*
	     * Region contains North Pole
	     */

	    return !lnx;
	}

    }
    return lnx;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrCreate --
 *
 *	This procedure creates and initializes an new geolinearray.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

GeoLnArr
GeoLnArrCreate(nLinesMax)
    unsigned nLinesMax;
{
    struct GeoLnArr *lnArrPtr;
    int n;

    lnArrPtr = (GeoLnArr)CKALLOC(sizeof(*lnArrPtr));
    lnArrPtr->descr = NULL;
    lnArrPtr->lines = NULL;
    GeoLnArrSetDescr(lnArrPtr, "");
    lnArrPtr->nLines = lnArrPtr->nLinesMax = 0;
    lnArrPtr->nPts = lnArrPtr->nMax = 0;
    lnArrPtr->lonMax = lnArrPtr->latMax = -INT_MAX;
    lnArrPtr->lonMin = lnArrPtr->latMin = INT_MAX;
    lnArrPtr->lines = NULL;
    if (nLinesMax == 0) {
	return lnArrPtr;
    }
    lnArrPtr->lines = (struct GeoLn **)CKALLOC(nLinesMax
	    * sizeof(struct GeoLn *));
    lnArrPtr->nLinesMax = nLinesMax;
    for (n = 0; n < nLinesMax; n++) {
	lnArrPtr->lines[n] = NULL;
    }
    return lnArrPtr;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrSetDescr --
 *
 *	This procedure sets the descriptor for a geolinearray.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnArrSetDescr(lnArrPtr, descr)
    struct GeoLnArr *lnArrPtr;
    CONST char *descr;
{
    lnArrPtr->descr = CKREALLOC(lnArrPtr->descr, strlen(descr) + 1);
    strcpy(lnArrPtr->descr, descr);
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrSetAlloc --
 *
 *	This procedure changes the number of geolines a geolinearray can store.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnArrSetAlloc(lnArrPtr, nLinesMax)
    struct GeoLnArr *lnArrPtr;
    unsigned nLinesMax;
{
    int n;		/* Loop index */

    if (lnArrPtr->nLinesMax == nLinesMax) {
	return;
    }
    for (n = nLinesMax; n < lnArrPtr->nLinesMax; n++) {
	/*
	 * Free excess lines
	 */

	GeoLnDestroy(lnArrPtr->lines[n]);
    }
    lnArrPtr->lines = (GeoLn *)CKREALLOC((char *)lnArrPtr->lines,
	    nLinesMax * sizeof(GeoLn));
    lnArrPtr->nLinesMax = nLinesMax;
    for (n = lnArrPtr->nLines; n < lnArrPtr->nLinesMax; n++) {
	/*
	 * Initialize additional lines.
	 */

	lnArrPtr->lines[n] = NULL;
    }

}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrAddLine --
 *
 *	This procedure copies a geoline onto the end of a geolinearray.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

int
GeoLnArrAddLine(ln, lnArrPtr)
    GeoLn ln;
    struct GeoLnArr *lnArrPtr;
{
    int nLines;			/* Initial number of lines in lnArrPtr */

    nLines = lnArrPtr->nLines;
    if (nLines + 1 > lnArrPtr->nLinesMax) {
	/*
	 * If lnArrPtr needs more space in its allocation, grow it by 25%
	 */

	GeoLnArrSetAlloc(lnArrPtr, ((lnArrPtr->nLinesMax + 4) * 5) / 4);
    }
    if ( !(lnArrPtr->lines[nLines] = GeoLnCreate(ln->nPts)) ) {
	return 0;
    }
    lnArrPtr->nPts += ln->nPts;
    lnArrPtr->nMax = (lnArrPtr->nMax > ln->nPts) ? lnArrPtr->nMax : ln->nPts;
    lnArrPtr->latMax
	= (lnArrPtr->latMax > ln->latMax) ? lnArrPtr->latMax : ln->latMax;
    lnArrPtr->lonMax
	= (lnArrPtr->lonMax > ln->lonMax) ? lnArrPtr->lonMax : ln->lonMax;
    lnArrPtr->latMin
	= (lnArrPtr->latMin < ln->latMin) ? lnArrPtr->latMin : ln->latMin;
    lnArrPtr->lonMin
	= (lnArrPtr->lonMin < ln->lonMin) ? lnArrPtr->lonMin : ln->lonMin;
    memcpy(lnArrPtr->lines[nLines]->pts, ln->pts, ln->nPts * sizeof(GeoPt));
    lnArrPtr->lines[nLines]->nPts = ln->nPts;
    lnArrPtr->lines[nLines]->lonMax = ln->lonMax;
    lnArrPtr->lines[nLines]->lonMin = ln->lonMin;
    lnArrPtr->lines[nLines]->latMax = ln->latMax;
    lnArrPtr->lines[nLines]->latMin = ln->latMin;
    lnArrPtr->nLines++;
    return 1;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrPutLine --
 *
 *	This procedure gives a geoline to a geolinearray.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

int
GeoLnArrPutLine(ln, lnArrPtr)
    GeoLn ln;
    struct GeoLnArr *lnArrPtr;
{
    int nLines;			/* Initial number of lines in lnArrPtr */

    nLines = lnArrPtr->nLines;
    if (nLines + 1 > lnArrPtr->nLinesMax) {
	/*
	 * If lnArrPtr needs more space in its allocation, grow it by 25%
	 */

	GeoLnArrSetAlloc(lnArrPtr, ((lnArrPtr->nLinesMax + 4) * 5) / 4);
    }
    lnArrPtr->nPts += ln->nPts;
    lnArrPtr->nMax = (lnArrPtr->nMax > ln->nPts) ? lnArrPtr->nMax : ln->nPts;
    lnArrPtr->latMax
	= (lnArrPtr->latMax > ln->latMax) ? lnArrPtr->latMax : ln->latMax;
    lnArrPtr->lonMax
	= (lnArrPtr->lonMax > ln->lonMax) ? lnArrPtr->lonMax : ln->lonMax;
    lnArrPtr->latMin
	= (lnArrPtr->latMin < ln->latMin) ? lnArrPtr->latMin : ln->latMin;
    lnArrPtr->lonMin
	= (lnArrPtr->lonMin < ln->lonMin) ? lnArrPtr->lonMin : ln->lonMin;
    lnArrPtr->lines[nLines] = ln;
    lnArrPtr->nLines++;
    return 1;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrGetDescr --
 *
 *	This is a member access function.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

char *
GeoLnArrGetDescr(lnArrPtr)
    struct GeoLnArr *lnArrPtr;
{
	return lnArrPtr->descr;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrGetLine --
 *
 *	This procedure retrieves a geoline from a geolinearray.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

GeoLn
GeoLnArrGetLine(lnArrPtr, n)
    struct GeoLnArr *lnArrPtr;
    unsigned n;
{
    return (n < lnArrPtr->nLines) ? lnArrPtr->lines[n] : NULL;
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrFree --
 *
 *	This procedure frees internal storage in a geolinearray.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnArrFree(lnArrPtr)
    struct GeoLnArr *lnArrPtr;
{
    int n;

    if ( !lnArrPtr ) {
	return;
    }
    for (n = 0; n < lnArrPtr->nLines; n++) {
	GeoLnDestroy(lnArrPtr->lines[n]);
    }
    CKFREE((char *)lnArrPtr->lines);
    CKFREE((char *)lnArrPtr->descr);
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrDestroy --
 *
 *	This procedure frees internal storage in a geolinearray and frees the
 *	geolinearray.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	See the user documentation.
 *
 *----------------------------------------------------------------------
 */

void
GeoLnArrDestroy(lnArrPtr)
    struct GeoLnArr *lnArrPtr;
{
    GeoLnArrFree(lnArrPtr);
    CKFREE((char *)lnArrPtr);
}

/*
 *----------------------------------------------------------------------
 *
 * GeoLnArrContainGeoPt --
 *
 *	This procedure determines whether any of the geolines in a geolinearray
 *	contain a geopoint.
 *
 * Results:
 *	See the user documentation.
 *
 * Side effects:
 *	None.
 *
 *----------------------------------------------------------------------
 */

int
GeoLnArrContainGeoPt(geoPt, lnArrPtr)
    GeoPt geoPt;
    struct GeoLnArr *lnArrPtr;
{
    int n;

    for (n = 0; n < lnArrPtr->nLines; n++) {
	if ( GeoLnContainGeoPt(geoPt, lnArrPtr->lines[n]) ) {
	    return 1;
	}
    }
    return 0;
}


syntax highlighted by Code2HTML, v. 0.9.1