[comp.lang.c] Need help with finding perimeters of bounded areas

miller@unicom.UUCP (Gregory S Miller) (10/12/89)

Consider an area A on an euclidean plane surface.  Inside A 
are an arbitrary number of lines that "fill out" A.  That is,
A is defined not by a perimeter, but rather as a filled object (
negative versus positive images...).

Now each line has a diameter D(i) where i is an subscript telling 
which line we`re referring.  It is perfectly reasonable that the lines
may overlap.  Since they overlap, one could use a lot of lines with
small diameters or a smaller number with larger diameters and still
produce the same area A.

Problem: Given the area A, find the perimeter - eliminate all lines
	 inside A and leave just the outline.  A is "given" by a
	 list of vectors with a diameter (ie. start/end point with
	 diameter). 

	 Just in case it was not already clear, the number, location
	 and diameter (D(i)) of each line is arbitrary so long as
	 the area A is filled out.

I do not know quite where to start on this problem.  I`ve never
run across "standard algorithms" which would be suited to this problem.

Ideas, text references, or the like would be helpful.
Thanks.

dhw@itivax.iti.org (David H. West) (10/13/89)

In article <135@unicom.UUCP> miller@unicom.UUCP (Gregory S Miller) writes:
>Consider an area A on an euclidean plane surface.  Inside A 
>are an arbitrary number of lines that "fill out" A.  That is,
>A is defined not by a perimeter, but rather as a filled object (
>negative versus positive images...).

>Problem: Given the area A, find the perimeter - eliminate all lines
>	 inside A and leave just the outline.  A is "given" by a
>	 list of vectors with a diameter (ie. start/end point with
>	 diameter). 

I imagine one could simplify one of the perimeter-marching convex hull
algorithms fairly straightforwardly to do this. See e.g.
Preparata and Shamos, Computational Geometry, Springer, 1985.

-David West        dhw@itivax.iti.org

decot@hpisod2.HP.COM (Dave Decot) (10/14/89)

One needs to know the width of the "lines", since otherwise an infinite
number of "lines" would be required (i.e., i is infinite), assuming that
A actually occupies non-zero area.

Dave

gwyn@smoke.BRL.MIL (Doug Gwyn) (10/14/89)

In article <135@unicom.UUCP> miller@unicom.UUCP (Gregory S Miller) writes:
>Consider an area A on an euclidean plane surface.  Inside A 
>are an arbitrary number of lines that "fill out" A.  That is,
>A is defined not by a perimeter, but rather as a filled object (
>negative versus positive images...).
>Now each line has a diameter D(i) where i is an subscript telling 
>which line we`re referring.  It is perfectly reasonable that the lines
>may overlap.  Since they overlap, one could use a lot of lines with
>small diameters or a smaller number with larger diameters and still
>produce the same area A.
>Problem: Given the area A, find the perimeter - eliminate all lines
>	 inside A and leave just the outline.  A is "given" by a
>	 list of vectors with a diameter (ie. start/end point with
>	 diameter). 
>	 Just in case it was not already clear, the number, location
>	 and diameter (D(i)) of each line is arbitrary so long as
>	 the area A is filled out.

If I understand the problem statement correctly, you can make a
preliminary processing pass wherein each given line is replaced
by its outline (rectangle of length = line's and width = D(i)),
then find the bounding polygon by the following algorithm that I
developed many years ago; it's not elegant but it works.

========== boundp.c: ==========
/*
	boundp -- bounding polygon of two-dimensional view

	last edit:	85/04/01	D A Gwyn

NOTES FOR MAINTAINER:

	This program is somewhat slow when operating on large input
	sets.  The following ideas should improve the speed, possibly
	at the expense of maximum data set size:

	(1)  Merge the Chop() and Build() phases.  Since Chop() examines
	segment endpoints anyway, it is in a good position to build the
	linked lists used by Search().

	(2)  Sort the input endpoints into a double tree on X and Y.
	This should cut Chop() from order N ^ 2 to N * log( N ).
*/
#ifndef	lint
static char	sccsid[] = "@(#)boundp.c	1.12";
#endif

#include	<math.h>
#include	<stdio.h>
#include	<string.h>
#ifdef	unix
#include	<varargs.h>
#else
#include	<stdarg.h>
#endif

#include	<std.h>

extern void	exit(), free();
extern pointer	malloc();
extern int	getopt();

#define datan2( y, x )	(atan2( (double)(y), (double)(x) ) * DEGRAD)

extern char	*optarg;

typedef struct
	{
	/* The following are float instead of double, to save space: */
	float		x;		/* X coordinate */
	float		y;		/* Y coordinate */
	}	coords; 		/* view-coordinates of point */

typedef struct segment
	{
	struct segment	*links; 	/* segment list thread */
	coords		sxy;		/* (X,Y) coordinates of start */
	coords		exy;		/* (X,Y) coordinates of end */
	}	segment;		/* entry in list of segments */

struct queue;				/* (for forward reference) */

typedef struct point
	{
	struct point	*linkp; 	/* point list thread */
	struct queue	*firstq;	/* NULL once output */
	coords		xy;		/* (X,Y) coordinates of point */
	}	point;			/* entry in list of points */

typedef struct queue
	{
	struct queue	*nextq; 	/* endpoint list thread */
	point		*endpoint;	/* -> endpt of ray from point */
	}	queue;			/* entry in list of endpoints */

static bool	Build(), Chop(), EndPoint(), GetArgs(), Input(), Mess(),
		Near(), Search(), Split(), Usage();
static coords	*Intersect();
static point	*LookUp(), *NewPoint(), *PutList();
static pointer	Alloc();
static queue	*Enqueue();
static void	Output(), Toss();

static bool	initial = true; 	/* false after first Output */
static bool	vflag = false;		/* set if "-v" option found */
static double	tolerance = 0.0;	/* point matching slop */
static double	tolsq = 0.0;		/* `tolerance' ^ 2 */
static point	*headp = NULL;		/* head of list of points */
static segment	seghead = { &seghead };	/* segment list head */


static bool
Usage() 				/* print usage message */
	{
	return
	Mess( "usage: boundp[ -i input][ -o output][ -t tolerance][ -v]"
	    );
	}


main( argc, argv )			/* "boundp" entry point */
	int	argc;			/* argument count */
	char	*argv[];		/* argument strings */
	{
	if ( !GetArgs( argc, argv ) )	/* process command arguments */
		return 1;

	if ( !Chop() )			/* input; chop into segments */
		return 2;

	if ( !Build() ) 		/* build linked lists */
		return 3;

	if ( !Search() )		/* output bounding polygon */
		return 4;

	return 0;			/* success! */
	}


static bool
GetArgs( argc, argv )			/* process command arguments */
	int		argc;		/* argument count */
	char		*argv[];	/* argument strings */
	{
	static bool	iflag = false;	/* set if "-i" option found */
	static bool	oflag = false;	/* set if "-o" option found */
	static bool	tflag = false;	/* set if "-t" option found */
	int		c;		/* option letter */

#ifdef	DEBUG
	(void)Mess( "\n\t\tGetArgs\n" );
#endif
	while ( (c = getopt( argc, argv, "i:o:t:v" )) != EOF )
		switch ( c )
			{
		case 'i':
			if ( iflag )
				return
				Mess( "too many -i options" );
			iflag = true;

			if ( strcmp( optarg, "-" ) != 0
			  && freopen( optarg, "r", stdin ) == NULL
			   )
				return
				Mess( "can't open \"%s\"", optarg );
			break;

		case 'o':
			if ( oflag )
				return
				Mess( "too many -o options" );
			oflag = true;

			if ( strcmp( optarg, "-" ) != 0
			  && freopen( optarg, "w", stdout ) == NULL
			   )
				return
				Mess( "can't create \"%s\"", optarg );
			break;

		case 't':
			if ( tflag )
				return
				Mess( "too many -t options" );
			tflag = true;

			if ( sscanf( optarg, "%le", &tolerance ) != 1 )
				return
				Mess( "bad tolerance: %s", optarg );
			tolsq = tolerance * tolerance;
			break;

		case 'v':
			vflag = true;
			break;

		case '?':
			return Usage(); /* print usage message */
			}

	return true;
	}


/*
	The following phase is required to chop up line segments that
	intersect; this guarantees that all polygon vertices will be
	found at segment endpoints.  Unfortunately, this is an N^2
	process.  It is also hard to do right; thus the elaborateness.
*/

static bool
Chop()					/* chop vectors into segments */
	{
	segment *inp;			/* -> input whole segment */

#ifdef	DEBUG
	(void)Mess( "\n\t\tChop\n" );
#endif
	/* Although we try to Alloc() when no input remains, it is more
	   efficient than repeatedly copying the input into *inp.     */
	while ( (inp = (segment *)Alloc( sizeof(segment) )) != NULL
	     && Input( inp )
	      ) {
		register segment	*segp;	/* segment list entry */
		segment 		piecehead;
						/* list of inp pieces */
		register segment	*pp;	/* -> pieces of `inp' */

		/* Reverse the segment if necessary to get endpoints
		   in left-to-right order; speeds up range check. */

		if ( inp->sxy.x > inp->exy.x )
			{
			coords	temp;		/* store for swapping */

#ifdef	DEBUG
			(void)Mess( "endpoints swapped" );
#endif
			temp = inp->sxy;
			inp->sxy = inp->exy;
			inp->exy = temp;
			}

		/* Split input and segment list entries into pieces. */

		inp->links = &piecehead;
		piecehead.links = inp;	/* initially, one piece */

		/* Note:  Although the split-off pieces of `segp' are
		   linked in front of `segp', that's okay since they
		   cannot intersect `inp' again!		      */
		for ( segp = seghead.links; segp != &seghead;
		      segp = segp->links
		    )
			/* Similarly, once a piece of `inp' intersects
			   `segp', we can stop looking at the pieces. */
			for ( pp = piecehead.links; pp != &piecehead;
			      pp = pp->links
			    )	{
				register coords *i;	/* intersectn */
	
				/* Be careful; `i' -> static storage
				   internal to Intersect().	      */
				i = Intersect( pp, segp );
				if ( i != NULL )
					{
					if ( !EndPoint( i, pp  )
					  && !Split( i, pp, &piecehead )
					  || !EndPoint( i, segp )
					  && !Split( i, segp, &seghead )
					   )	/* out of memory */
						return false;

					break;	/* next `segp' */
					}
				}
#ifdef	DEBUG
		(void)Mess( "new input pieces:" );
		for ( pp = piecehead.links; pp != &piecehead;
		      pp = pp->links
		    )
			(void)Mess( "\t(%g,%g)->(%g,%g)",
				    (double)pp->sxy.x,
				    (double)pp->sxy.y,
				    (double)pp->exy.x,
				    (double)pp->exy.y
				  );
		(void)Mess( "other segments:" );
		for ( segp = seghead.links; segp != &seghead;
		      segp = segp->links
		    )
			(void)Mess( "\t(%g,%g)->(%g,%g)",
				    (double)segp->sxy.x,
				    (double)segp->sxy.y,
				    (double)segp->exy.x,
				    (double)segp->exy.y
				  );
#endif

		/* Put input segment pieces into segment list. */

		/* Note:  It is better to scan the pieces to find the
		   tail link than to waste storage on double links. */
		for ( pp = piecehead.links; pp->links != &piecehead;
		      pp = pp->links
		    )
			;
		pp->links = seghead.links;	/* last-piece link */
		seghead.links = piecehead.links;	/* add pieces */
		}
	Toss( (pointer)inp );		/* unused segment at EOF */

	return inp != NULL;
	}


static bool
Split( p, oldp, listh ) 		/* split segment in two */
	coords			*p;	/* -> break point */
	register segment	*oldp;	/* -> segment to be split */
	register segment	*listh; /* -> list to attach tail to */
	{
	register segment	*newp;	/* -> new list entry */

#ifdef	DEBUG
	(void)Mess( "split (%g,%g)->(%g,%g) at (%g,%g)",
		    (double)oldp->sxy.x, (double)oldp->sxy.y,
		    (double)oldp->exy.x, (double)oldp->exy.y,
		    (double)p->x, (double)p->y
		  );
#endif
	if ( (newp = (segment *)Alloc( sizeof(segment) )) == NULL )
		return false;		/* out of heap space */
	newp->links = listh->links;
	newp->sxy = *p;
	newp->exy = oldp->exy;
	oldp->exy = *p; 		/* old entry, new endpoint */
	listh->links = newp;		/* attach to list head */

	return true;
	}


static bool
Build() 				/* build linked lists */
	{
	register segment	*listp; /* -> segment list entry */
	segment 		*deadp; /* -> list entry to be freed */

#ifdef	DEBUG
	(void)Mess( "\n\t\tBuild\n" );
#endif
	/* When we are finished, `seghead' will become invalid. */
	for ( listp = seghead.links; listp != &seghead;
	      deadp = listp, listp = listp->links, Toss( (pointer)deadp )
	    )	{
		register point	*startp, *endp; /* -> segment endpts */

		if ( (startp = PutList( &listp->sxy )) == NULL
		  || (endp   = PutList( &listp->exy )) == NULL
		  || Enqueue( startp, endp ) == NULL
		  || Enqueue( endp, startp ) == NULL
		   )
			return false;	/* out of heap space */
		}

	return	true;
	}


static bool
Search()				/* output bounding polygon */
	{
	double		from;		/* backward edge direction */
	register point	*currentp;	/* -> current (start) point */
	point		*previousp;	/* -> previous vertex point */

#ifdef	DEBUG
	(void)Mess( "\n\t\tSearch\n" );
#endif
	if ( headp == NULL )
		return true;		/* trivial case */

	/* Locate the lowest point; this is the first polygon vertex. */
	{
	float		miny;		/* smallest Y coordinate */
	register point	*listp; 	/* -> next point in list */

	currentp = listp = headp;
	miny = currentp->xy.y;
	while ( (listp = listp->linkp) != NULL )
		if ( listp->xy.y < miny )
			{
			miny = listp->xy.y;
			currentp = listp;
			}
	}
	previousp = NULL;		/* differs from currentp! */
	from = 270.0;

	/* Output point and look for next CCW point on periphery. */

	for ( ; ; )
		{
		coords		first;	/* first output if -v */
		double		mindir; /* smallest from->to angle */
		register point	*nextp; /* -> next perimeter point */
		queue		*endq;	/* -> endpoint queue entry */

#ifdef	DEBUG
		(void)Mess( "from %g", from );
#endif
		endq = currentp->firstq;
		if ( endq == NULL )
			{
			if ( vflag && !initial )
				Output( &first );	/* closure */

			return true;	/* been here before => done */
			}

		Output( &currentp->xy );	/* found vertex */
		if ( vflag && initial )
			{
			initial = false;
			first = currentp->xy;	/* save for closure */
			}

		/* Find the rightmost forward edge endpoint. */

		mindir = 362.0;
		do	{
			double		to;	/* forward edge dir */
			double		diff;	/* angle from->to */
			register point	*endp;	/* -> endpoint */

			endp = endq->endpoint;
			if ( endp == previousp	/* don't double back! */
			  || endp == currentp	/* don't stay here! */
			   )
				continue;

			/* Note: it would be possible to save some calls
			   to atan2 by being clever about quadrants.  */
			to = datan2( endp->xy.y - currentp->xy.y,
				     endp->xy.x - currentp->xy.x
				   );
#ifdef	DEBUG
			(void)Mess( "to %g", to );
#endif
			diff = to - from;
			/* Note: Exact 360 (0) case is not supposed to
			   happen, but this algorithm copes with it.  */
			while ( diff <= 0.0 )
				diff += 360.0;

			if ( diff < mindir )
				{
#ifdef	DEBUG
				(void)Mess( "new min %g", diff );
#endif
				mindir = diff;
				nextp = endp;
				}
			} while ( (endq = endq->nextq) != NULL );

		if ( mindir > 361.0 )
			return Mess( "degenerate input" );

		currentp->firstq = NULL;	/* "visited" */
		previousp = currentp;
		currentp = nextp;

		from += mindir + 180.0; /* reverse of saved "to" */
		/* The following is needed only to improve accuracy: */
		while ( from > 360.0 )
			from -= 360.0;	/* reduce to range [0,360) */
		}
	}


static point *
PutList( coop ) 			/* return -> point in list */
	register coords *coop;		/* -> coordinates */
	{
	register point	*p;		/* -> list entry */

	p = LookUp( coop );		/* may already be there */
	if ( p == NULL )		/* not yet in list */
		{			/* start new point group */
#ifdef	DEBUG
		(void)Mess( "new point group (%g,%g)",
			    (double)coop->x, (double)coop->y
			  );
#endif
		p = NewPoint( coop );
/*		if ( p == NULL )
/*			return NULL;	/* out of heap space */
		}

	return p;			/* -> point list entry */
	}


static point *
LookUp( coop )				/* find point group in list */
	register coords *coop;		/* -> coordinates */
	{
	register point	*p;		/* -> list members */

	for ( p = headp; p != NULL; p = p->linkp )
		if ( Near( coop, &p->xy ) )
			{
#ifdef	DEBUG
			(void)Mess( "found (%g,%g) in list",
				    (double)coop->x, (double)coop->y
				  );
#endif
			return p;	/* found a match */
			}

	return NULL;			/* not yet in list */
	}


static point *
NewPoint( coop )			/* add point to list */
	register coords *coop;		/* -> coordinates */
	{
	register point	*newp;		/* newly allocated point */

	newp = (point *)Alloc( sizeof(point) );
	if ( newp == NULL )
		return NULL;

#ifdef	DEBUG
	(void)Mess( "add point (%g,%g)",
		    (double)coop->x, (double)coop->y
		  );
#endif
	newp->linkp = headp;
	newp->firstq = NULL;		/* empty endpoint queue */
	newp->xy = *coop;		/* coordinates */
	return headp = newp;
	}


static queue *
Enqueue( addp, startp ) 		/* add to endpoint queue */
	register point	*addp;		/* -> point being queued */
	register point	*startp;	/* -> point owning queue */
	{
	register queue	*newq;		/* new queue element */

	newq = (queue *)Alloc( sizeof(queue) );
	if ( newq == NULL )
		return NULL;

#ifdef	DEBUG
	(void)Mess( "enqueue (%g,%g) on (%g,%g)",
		    (double)addp->xy.x, (double)addp->xy.y,
		    (double)startp->xy.x, (double)startp->xy.y
		  );
#endif
	newq->nextq = startp->firstq;
	newq->endpoint = addp;
	return startp->firstq = newq;
	}


static coords *
Intersect( a, b )			/* determine intersection */
	register segment	*a, *b; /* segments being tested */
	{
	double			det;	/* determinant, 0 if parallel */
	double			xaeas, xbebs, yaeas, ybebs;
					/* coordinate differences */

	/* First perform range check, to eliminate most cases.	Note
	   that segments point left-to-right even after splitting. */

	if ( a->sxy.x > 		/* a left */
	     b->exy.x + tolerance	/* b right */
	  || a->exy.x < 		/* a right */
	     b->sxy.x - tolerance	/* b left */
	  || Min( a->sxy.y, a->exy.y ) >		/* a bottom */
	     Max( b->sxy.y, b->exy.y ) + tolerance	/* b top */
	  || Max( a->sxy.y, a->exy.y ) <		/* a top */
	     Min( b->sxy.y, b->exy.y ) - tolerance	/* b bottom */
	   )	{
#ifdef	DEBUG
		(void)Mess( "ranges don't intersect" );
#endif
		return NULL;		/* can't intersect */
		}

	/* Passed quick check, now comes the hard part. */

	xaeas = (double)a->exy.x - (double)a->sxy.x;
	xbebs = (double)b->exy.x - (double)b->sxy.x;
	yaeas = (double)a->exy.y - (double)a->sxy.y;
	ybebs = (double)b->exy.y - (double)b->sxy.y;

	det = xbebs * yaeas - xaeas * ybebs;

	{
	double	norm;			/* norm of coefficient matrix */
	double	t;			/* test value for norm */

	norm = 0.0;
	if ( (t = Abs( xaeas )) > norm )
		norm = t;
	if ( (t = Abs( xbebs )) > norm )
		norm = t;
	if ( (t = Abs( yaeas )) > norm )
		norm = t;
	if ( (t = Abs( ybebs )) > norm )
		norm = t;

#define EPSILON 1.0e-06 		/* relative `det' size thresh */
	if ( Abs( det ) <= EPSILON * norm * norm )
		{
#ifdef	DEBUG
		(void)Mess( "parallel: det=%g, norm=%g", det, norm );
#endif
		return NULL;		/* parallels don't intersect */
		}
#undef	EPSILON
	}
	{
	/* `p' must be static; Intersect returns a pointer to it! */
	static coords	p;		/* point of intersection */
	double		lambda, mu;	/* segment parameters */
	double		onemmu; 	/* 1.0 - mu, for efficiency */
	double		xbsas, ybsas;	/* more coord differences */

	xbsas = (double)b->sxy.x - (double)a->sxy.x;
	ybsas = (double)b->sxy.y - (double)a->sxy.y;

	mu = (xbebs * ybsas - xbsas * ybebs) / det;
	onemmu = 1.0 - mu;
	p.x = onemmu * a->sxy.x + mu * a->exy.x;
	p.y = onemmu * a->sxy.y + mu * a->exy.y;
	if ( (onemmu < 0.0 || mu < 0.0) && !EndPoint( &p, a ) )
		{
#ifdef	DEBUG
		(void)Mess( "intersect off (%g,%g)->(%g,%g): mu=%g",
			    (double)a->sxy.x, (double)a->sxy.y,
			    (double)a->exy.x, (double)a->exy.y,
			    mu
			  );
#endif
		return NULL;		/* not in segment *a */
		}

	lambda = (xaeas * ybsas - xbsas * yaeas) / det;
	if ( (lambda > 1.0 || lambda < 0.0) && !EndPoint( &p, b ) )
		{
#ifdef	DEBUG
		(void)Mess( "intersect off (%g,%g)->(%g,%g): lambda=%g",
			    (double)b->sxy.x, (double)b->sxy.y,
			    (double)b->exy.x, (double)b->exy.y,
			    lambda
			  );
#endif
		return NULL;		/* not in segment *b */
		}

#ifdef	DEBUG
	(void)Mess( "intersection is (%g,%g): mu=%g lambda=%g",
		    (double)p.x, (double)p.y, mu, lambda
		  );
#endif
	return &p;
	}
	}


static bool
EndPoint( p, segp )			/* check for segment endpoint */
	register coords 	*p;	/* -> point being tested */
	register segment	*segp;	/* -> segment */
	{
#ifdef	DEBUG
	if ( Near( p, &segp->sxy ) || Near( p, &segp->exy ) )
		(void)Mess( "(%g,%g) is endpt of (%g,%g)->(%g,%g)",
			    (double)p->x, (double)p->y,
			    (double)segp->sxy.x, (double)segp->sxy.y,
			    (double)segp->exy.x, (double)segp->exy.y
			  );
#endif
	return Near( p, &segp->sxy ) || Near( p, &segp->exy );
	}


static bool
Near( ap, bp )				/* check if within tolerance */
	register coords *ap, *bp;	/* -> points being checked */
	{
	double		xsq, ysq;	/* dist between coords ^ 2 */

	/* Originally this was an abs value test; this is neater. */

	xsq = ap->x - bp->x;
	xsq *= xsq;
	ysq = ap->y - bp->y;
	ysq *= ysq;

#ifdef	DEBUG
	if ( xsq + ysq <= tolsq )
		(void)Mess( "(%g,%g) is near (%g,%g)",
			    (double)ap->x, (double)ap->y,
			    (double)bp->x, (double)bp->y
			  );
#endif
	return xsq + ysq <= tolsq;
	}


static pointer
Alloc( size )				/* allocate storage from heap */
	unsigned		size;	/* # bytes required */
	{
	register pointer	ptr;	/* -> allocated storage */

	if ( (ptr = malloc( size )) == NULL )
		(void)Mess( "out of memory" );

	return ptr;			/* (may be NULL) */
	}


static void
Toss( ptr )				/* return storage to heap */
	register pointer	ptr;	/* -> allocated storage */
	{
	if ( ptr != NULL )
		free( ptr );
	}


/*VARARGS0*/
static bool
#ifdef	unix
Mess( va_alist )			/* print error message */
	va_dcl				/* format, optional arguments */
#else	/* ANSI C */
Mess( fmt, ... )			/* print error message */
	register char	*fmt;		/* format */
#endif
	{
	va_list		ap;		/* for accessing arguments */
#ifdef	unix
	register char	*fmt;		/* format */

	va_start( ap );
	fmt = va_arg( ap, char * );
#else
	va_start( ap, fmt );
#endif
	(void)fflush( stdout );
	(void)fputs( "boundp: ", stderr );
	(void)vfprintf( stderr, fmt, ap );
	(void)fputc( '\n', stderr );

	va_end( ap );

	return false;
	}


static bool
Input( inp )				/* input stroke record */
	register segment	*inp;	/* -> input segment */
	{
	char			inbuf[82];	/* record buffer */

	while ( fgets( inbuf, (int)sizeof inbuf, stdin ) != NULL )
		{			/* scan input record */
		register int	cvt;	/* # fields converted */

#ifdef	DEBUG
		(void)Mess( "input: %s", inbuf );
#endif
		cvt = sscanf( inbuf, " %e %e %e %e",
			      &inp->sxy.x, &inp->sxy.y,
			      &inp->exy.x, &inp->exy.y
			    );

		if ( cvt == 0 )
			continue;	/* skip color, comment, etc. */

		if ( cvt == 4 )
			return true;	/* successfully converted */

		(void)Mess( "bad input: %s", inbuf );
		exit( 5 );		/* return false insufficient */
		}

	return false;			/* EOF */
	}


static void
Output( coop )				/* dump polygon vertex coords */
	register coords *coop;		/* -> coords to be output */
	{
	static coords	last;		/* previous *coop */

	if ( vflag )
		{
		if ( !initial )
			(void)printf( "%g %g %g %g\n",
				      (double)last.x, (double)last.y,
				      (double)coop->x, (double)coop->y
				    );

		last = *coop;		/* save for next start point */
		}
	else
		(void)printf( "%g %g\n",
			      (double)coop->x, (double)coop->y
			    );
#ifdef	DEBUG
	(void)Mess( "output: %g %g",
		    (double)coop->x, (double)coop->y
		  );
#endif
	}
========== std.h: ==========
/*
	std.h -- Douglas A. Gwyn's standard C programming definitions

	Prerequisite:	<math.h> (if you invoke Round())

	last edit:	89/09/28	D A Gwyn

	SCCS ID:	@(#)std.h	1.32

	The master source file is to be modified only by Douglas A. Gwyn
	<Gwyn@BRL.MIL>.  When installing a VLD/VMB software distribution,
	this file may need to be tailored slightly to fit the target system.
	Usually this just involves enabling some of the "kludges for deficient
	C implementations" at the end of this file.
*/

#ifndef	VLD_STD_H_
#define	VLD_STD_H_			/* once-only latch */

/* Extended data types */

typedef int	bool;			/* Boolean data */
#define 	false	0
#define 	true	1

typedef int	bs_type;		/* 3-way "bug/status" result type */
#define 	bs_good	1
#define 	bs_bad	0
#define 	bs_ugly	(-1)

/* ANSI C definitions */

/* Defense against some silly systems defining __STDC__ to random things. */
#ifdef STD_C
#undef STD_C
#endif
#ifdef __STDC__
#if __STDC__ > 0
#define	STD_C	__STDC__		/* use this instead of __STDC__ */
#endif
#endif

#ifdef STD_C
typedef void	*pointer;		/* generic pointer */
#else
typedef char	*pointer;		/* generic pointer */
#define	const		/* nothing */	/* ANSI C type qualifier */
/* There really is no substitute for the following, but these might work: */
#define	signed		/* nothing */	/* ANSI C type specifier */
#define	volatile	/* nothing */	/* ANSI C type qualifier */
#endif

#ifndef EXIT_SUCCESS
#define	EXIT_SUCCESS	0
#endif

#ifndef EXIT_FAILURE
#define	EXIT_FAILURE	1
#endif

#ifndef NULL
#define NULL	0			/* null pointer constant, all types */
#endif

/* Universal constants */

#define DEGRAD	57.2957795130823208767981548141051703324054724665642
					/* degrees per radian */
#define	E	2.71828182845904523536028747135266249775724709369996
					/* base of natural logs */
#define	GAMMA	0.57721566490153286061
					/* Euler's constant */
#define LOG10E	0.43429448190325182765112891891660508229439700580367
					/* log of e to the base 10 */
#define PHI	1.618033988749894848204586834365638117720309180
					/* golden ratio */
#define PI	3.14159265358979323846264338327950288419716939937511
					/* ratio of circumf. to diam. */

/* Useful macros */

/*
	The comment "UNSAFE" means that the macro argument(s) may be evaluated
	more than once, so the programmer must realize that the macro doesn't
	quite act like a genuine function.  This matters only when evaluating
	an argument produces "side effects".
*/

/* arbitrary numerical arguments and value: */
#define Abs( x )	((x) < 0 ? -(x) : (x))			/* UNSAFE */
#define Max( a, b )	((a) > (b) ? (a) : (b))			/* UNSAFE */
#define Min( a, b )	((a) < (b) ? (a) : (b))			/* UNSAFE */

/* floating-point arguments and value: */
#define Round( d )	floor( (d) + 0.5 )		/* requires <math.h> */

/* arbitrary numerical arguments, integer value: */
#define	Sgn( x )	((x) == 0 ? 0 : (x) > 0 ? 1 : -1)	/* UNSAFE */

/* string arguments, boolean value: */
#ifdef gould	/* UTX-32 2.0 compiler has problems with "..."[] */
#define	StrEq( a, b )	(strcmp( a, b ) == 0)			/* UNSAFE */
#else
#define	StrEq( a, b )	(*(a) == *(b) && strcmp( a, b ) == 0)	/* UNSAFE */
#endif

/* array argument, integer value: */
#define	Elements( a )	(sizeof a / sizeof a[0])

/* integer (or character) arguments and value: */
#define fromhostc( c )	(c)		/* map host char set to ASCII */
#define tohostc( c )	(c)		/* map ASCII to host char set */
#define tonumber( c )	((c) - '0')	/* convt digit char to number */
#define todigit( n )	((n) + '0')	/* convt digit number to char */

/* Kludges for deficient C implementations */

/*#define	strchr	index		/* 7th Edition UNIX, 4.2BSD */
/*#define	strrchr	rindex		/* 7th Edition UNIX, 4.2BSD */
/*#define	void	int		/* K&R Appendix A followers */

#endif	/* VLD_STD_H_ */

Howard.Spindel@p8.f14.n105.z1.FIDONET.ORG (Howard Spindel) (10/15/89)

> From: miller@unicom.UUCP (Gregory S Miller)
> Date: 12 Oct 89 02:37:36 GMT
> Organization: Science Computer Center, College of Marin,
> Kentfield CA
> Message-ID: <135@unicom.UUCP>
> Newsgroups: comp.lang.c,comp.lang.lisp,comp.ai
>
>
>
> Consider an area A on an euclidean plane surface.  Inside A
> are an arbitrary number of lines that "fill out" A.  That
> is,
> A is defined not by a perimeter, but rather as a filled
> object (
> negative versus positive images...).
>
> Now each line has a diameter D(i) where i is an subscript
> telling
> which line we`re referring.  It is perfectly reasonable that
> the lines
> may overlap.  Since they overlap, one could use a lot of
> lines with
> small diameters or a smaller number with larger diameters
> and still
> produce the same area A.
>
> Problem: Given the area A, find the perimeter - eliminate
> all lines
> inside A and leave just the outline.  A is "given"
> by a
> list of vectors with a diameter (ie. start/end
> point with
> diameter).
>
You don't say whether the area has any concavities or not.  If not,
the following brute force idea may work - I didn't spend a whole lot
of time thinking it through.
Make 4 lists.  List 1 is the minimum X coordinate seen for a given Y,
List 2 is the maximum X coordinate for a given Y, List 3 is the
minimum Y coordinate seen for a given X, List 4 is the maximum Y
coordinate seen for a given X.  Your arrays will have to be large
enough to hold an entry for each possible X and Y intercept.
Scan convert each line.  For each point encountered while scan
converting compare against the minima/maxima saved in the tables
and replace if new minima/maxima are found.
When all lines have been scan converted the lists should contain
the coordinates of each point along the edge of the desired area.
 
Eschew Sesquipedalian Obfuscation!
Usenet:     Howard.Spindel@p8.f14.n105.z1.FIDONET.ORG
Fidonet:    1:105/14.8


--  
Howard Spindel - via FidoNet node 1:105/14
	    UUCP: ...!{uunet!oresoft, tektronix!reed}!busker!14.8!Howard.Spindel
	    ARPA: Howard.Spindel@p8.f14.n105.z1.FIDONET.ORG