Computing Common Tangents Without a Separating Line

This is the source code that accompanies the 1995 WADS paper, "Computing common tangents without a separating line," by David Kirkpatrick and Jack Snoeyink ( nosep.ps.gz). It is written for the SGI, so it uses the GL graphics library. (It was also written to fit on two 8.5x11" pages, so is not a model of programming style.)

I need to find the time to put it into a better format.

main.h

/* main.h   Jack Snoeyink    Common tangent without a separating line
 *
 */

#pragma once

#include 
#include 
#include 

#define MAXPTS 2000              /* Maximum number of points per polyline */
#define EPSILON 1.0e-13         /* Approximation of zero */
#define SQRhf M_SQRT1_2
#define TWO_PI 2*M_PI

typedef double COORD;
typedef COORD POINT[2];         /* Most data is Cartesian points */
typedef COORD HOMOG[3];         /* Some partial calculations are Homogeneous */

#define XX 0
#define YY 1
#define WW 2

typedef struct Polygon { 
  int n,                        /* Number of vertices in polygon */
  st, end,			/* Tangent is in $(\st, \end]$ */
  tent,                         /* Index of tentative refinement if $\tent \ne \st$ */
  ccw,				/* 1 = ccw  -1 = cw */
  wrap;                         /* Boolean indicates wraparound */
  HOMOG tang;			/* Tangent $\tau_\tent$ when refined ($\tent\ne\st$) */
  POINT v[MAXPTS];
} Polygon;

#define DET2x2(p, q, i, j) ((p)[i]*(q)[j] - (p)[j]*(q)[i]) /* Determinants */
#define DET2(p, q) DET2x2(p,q, XX,YY) 
#define DET3C(p, q, r) DET2(q,r) - DET2(p,r) + DET2(p,q)

#define CROSSPROD_2CCH(p, q, r) /* 2-d Cartesian to Homogeneous cross product */\
  (r)[XX] = - (q)[YY] + (p)[YY];\
  (r)[YY] =   (q)[XX] - (p)[XX]; (r)[WW] = DET2(p,q); 

#define CROSSPROD_2SCCH(s, p, q, r) /* 2-d cartesian to homog cross product with sign */\
  (r)[XX] = s * (- (q)[YY] + (p)[YY]); \
  (r)[YY] = s * (  (q)[XX] - (p)[XX]);  (r)[WW] = s * DET2(p,q);

#define DOTPROD_2CH(p, q)       /* 2-d Cartesian to Homogeneous dot product */\
  ((q)[WW] + (p)[XX]*(q)[XX] + (p)[YY]*(q)[YY])

#define ASSIGN_H(p, op, q)	/* Homogeneous assignment */\
  (p)[WW] op (q)[WW];  (p)[XX] op (q)[XX];  (p)[YY] op (q)[YY];

#define LEFT(x) (pr(x, " LEFT=%g ") > EPSILON) /* Sidedness tests */
#define RIGHT(x) (pr(x, " RT=%g ") < -EPSILON)
#define LEFT_PL(p, l) LEFT(DOTPROD_2CH(p, l))
#define RIGHT_PL(p, l) RIGHT(DOTPROD_2CH(p, l))
#define LEFT_PPP(p, q, r) LEFT(DET3C(p, q, r))
#define RIGHT_PPP(p, q, r) RIGHT(DET3C(p, q, r))

extern int gr_verbose, verbose; /* debugging */
extern double pr();

nosep.c

/* nosep.c   Jack Snoeyink     Common tangent without a separating line
 *
 */
#include "main.h"

#define Pv(m) P->v[m]           /* Indexing into polygon vertices */
#define Qv(m) Q->v[m]
#define Pvr(m) P->v[(m) % P->n] /* Indexing mod $n$ */
#define Qvr(m) Q->v[(m) % Q->n]

#define CCW(x) (x->ccw == 1) /* Is $x$ oriented counterclockwise? */
#define DONE(x) ((x->end - x->tent) == x->ccw) /* Any candidates left? */
#define REFINED(x) (x->st != x->tent) /* Is $x$ refined? */

#define NO_DISC     -1           /* Actions in \op Refine() */
#define DISC_END     0
#define DISC_START   1


void Refine(P, Q)
     Polygon *P, *Q;            /* We refine polygon $P$ checking against $Q$.  
We can assume that more than one candidates exist in $(\tentP,\enP\,]$
and the invariants hold. */
{
  HOMOG q0pm, mtang;
  register int mid, left_base, action = NO_DISC;
  register COORD *pm, *pm1, *qend, *qt;
  
  mid = P->tent + (P->end - P->tent) / 2; /* Check \mid\ point.
                                            Round towards \tentP\ */
  pm = Pvr(mid);
  pm1 = Pvr(mid + P->ccw);
  CROSSPROD_2SCCH(P->ccw, pm, pm1, mtang); /* Generate $\tau_\mid$ */
  CROSSPROD_2CCH(Qv(0), pm, q0pm);
  left_base = RIGHT_PL(Pv(0), q0pm); 


  /* debugging */
#define RADIUS 0.02
  if (verbose) {
    printf("Refine %x %x %3d %3d %3d %3d \n", 
	   P, Q, P->st, P->tent, mid, P->end);
    fflush(stdout);
  }
  if (gr_verbose) {
    color(YELLOW);
    drawLine(mtang);
    color(GREEN); circle(Pvr(P->st), 1.0*RADIUS);
    color(BLUE); circle(pm, 1.5*RADIUS);

    color(RED); circle(Pvr(P->end), 2.0*RADIUS);
    displayAll();    } /**/
  
  if (REFINED(Q) && !LEFT_PL(pm, Q->tang)) {
    qt = Qvr(Q->tent);
    if (CCW(Q) ^ LEFT_PPP(Pv(0), qt, pm)) /* Check $\sigma_\tentQ$ */
      Q->st = Q->tent;          /* Certify tentative to $Q$ */
    else {
      Q->end = Q->tent; 
      Q->tent = Q->st;          /* Revoke tentative to $Q$ */
      P->st = P->tent;	/* Certify tentatve on $P$ (if refined) */
    }
  }

  qend = Qvr(Q->end);
  qt = Qvr(Q->tent);

  if (P->wrap && (left_base ^ (mid > P->n))) /* Is $P$ wrapped around? */
    action = !left_base;
  else if (!LEFT_PL(Qv(0), mtang)) /* Can we be tangent w.r.t $q_0$? */
    action = left_base;
  else if (!LEFT_PL(qend, mtang)) /* Be tangent w.r.t $q_\enQ$? */
    action = LEFT_PL(qend, q0pm);
  else if (REFINED(Q) && !LEFT_PL(qt, mtang)) /* Be tangent w.r.t $q_\tentQ$? */
    action = LEFT_PL(qt, q0pm);
  
  if (action == NO_DISC) {
    P->tent = mid;              /* We tentatively refine at \mid */        
    ASSIGN_H(P->tang, =, mtang)
    }
  else if (CCW(P) ^ action)
    P->st = P->tent = mid;      /* A discard at \stP\ occurred */
  else
    P->end = mid;               /* A discard at \enP\ occurred */
}




void Tang(P, Q)
     Polygon *P, *Q;
{                               /* Compute a tangent from $P$ to $Q$ */
  register int n1 = Q->n - 1;

  P->ccw = 1; P->st = P->tent = 0; P->end = P->n; /* Initialize $P$ */
  CROSSPROD_2CCH(Pv(0), Pv(1), P->tang);
  if (P->wrap = LEFT_PL(Qv(0), P->tang)) /* Wrap $P$ initially */
    { P->tent = P->n; P->end += P->n; }

  Q->ccw = -1; Q->st = Q->tent = Q->n; Q->end = 0; /* Initialize $Q$ */
  CROSSPROD_2CCH(Qv(n1), Qv(0), Q->tang);
  if (Q->wrap = LEFT_PL(Pv(0), Q->tang)) /* Wrap $Q$ initially */
    Q->st += Q->n; 

  if (gr_verbose) 
    displayAll();/**/

  while (!DONE(P) || !DONE(Q)) {
    if (!DONE(P)) Refine(P, Q);
    if (!DONE(Q)) Refine(Q, P);
  }                             /* Finished. $\enQ$ and $\enP$ indicate tangent */
}

main.c

/* main.c Jack Snoeyink  
 * tangent without a separating line
 */

#include "main.h"
#include 

#define wLEFT    -1.0
#define wRIGHT    3.1
#define wBOTTOM  -1.0
#define wTOP      1.0
#define OFFSET   0.1
#define RADIUS 0.02

#define Pv(m) P->v[m]		/* Indexing into polygon vertices */
#define Qv(m) Q->v[m % Q->n]
#define Pvr(m) P.v[(m) % P.n]	/* Indexing vertices mod P.n */
#define Qvr(m) Q.v[(m) % Q.n]
  
#define LINCOMB_2C(a, p, b, q, r) /* 2-d cartesian linear combination */\
  (r)[XX] = (a) * (p)[XX] + (b) * (q)[XX];\
  (r)[YY] = (a) * (p)[YY] + (b) * (q)[YY];

extern void Tang();

Polygon P, Q;	/* convex polygons */

float offp, offq;  /* angle offsets for generation */
int verbose, gr_verbose;
int np, nq;
int randomize;
float rad1, rad2, rad3, rad4;

float rsize[8] = {1.05, 1.0, .95, .9, .8, .5, .1, .05};

void fatalError(msg, var)
     char *msg, *var;
{
  fprintf(stderr, msg, var);
  fflush(stderr);
  exit(1);
}

double pr(x, s)
     double x;
     char *s;
{
  if (verbose)
    {
      printf(s, x);
      fflush(stdout);
    }
  return (x);
}

void circle(p, r)
     POINT p;
     double r;
{
  int i;
  POINT tmp;
  static POINT off[8] = {1,0, SQRhf,SQRhf, 0,1, -SQRhf,SQRhf, 
			   -1,0, -SQRhf,-SQRhf, 0,-1, SQRhf,-SQRhf};
  bgnclosedline();
  for (i = 0; i < 8; i++)
    {
      LINCOMB_2C(1.0, p, r, off[i], tmp);
      v2d(tmp);
    }
  endclosedline();
}

void drawSeg(p, q)
     POINT p, q;
{
  bgnline();
  v2d(p);
  v2d(q);
  endline();
  gflush();
}


void drawLine(l)
     HOMOG l;
{
  POINT p, q;
  
  if (fabs(l[XX]) > fabs(l[YY]))
    {
      p[XX] = (-l[WW]-wBOTTOM*l[YY]) / l[XX];
      p[YY] = wBOTTOM;
      q[XX] = (-l[WW]-wTOP*l[YY]) / l[XX];
      q[YY] = wTOP;
    }
  else
    {
      p[XX] = wLEFT;
      p[YY] = (-l[WW]-wLEFT*l[XX]) / l[YY];
      q[XX] = wRIGHT;
      q[YY] = (-l[WW]-wRIGHT*l[XX]) / l[YY];
    }
  drawSeg(p, q);
}

void drawPoly(Q)
     Polygon *Q;
{
  register int i;
  
  if (verbose)
    printf(" \n", 
	   Q, Q->ccw, Q->st, Q->tent, Q->end);
  
  if (gr_verbose)
    {
      bgnclosedline();
      for (i = 0; i < Q->n; i++)
	v2d(Qv(i)); 
      endclosedline();
      fflush(stdout);
      color(BLUE); circle(Qv(Q->tent), 1.5*RADIUS);
      color(GREEN); circle(Qv(Q->st), 1*RADIUS);
      color(RED); circle(Qv(Q->end), 2.0*RADIUS);
    }
}

#define REFINED(x) (x.st != x.tent) /* is x refined? */

void displayAll()
{
  if (gr_verbose)
    {
      color(BLUE);
      if (REFINED(P)) drawLine(P.tang);
      if (REFINED(Q)) drawLine(Q.tang);
      color(MAGENTA);
    }
  drawPoly(&P);
  
  if (gr_verbose)
    color(CYAN);
  drawPoly(&Q);
  if (gr_verbose)
    {
      swapbuffers();
      color(BLACK);
      clear();
      if (verbose>1) 
	while (getchar() <= 0);
    }
}

void drawInit()
{
  foreground();
  prefposition(600,1000,200,400);
  winopen("nosep");
  wintitle("Tangents without a separating line");
  mmode(MVIEWING);
  doublebuffer();
  ortho2(wLEFT-OFFSET, wRIGHT+OFFSET, wBOTTOM-OFFSET, wTOP+OFFSET);
  gconfig();
  color(BLACK);
  clear();
  swapbuffers();
  clear();
}

void getPoints(P, n, ctrx, ctry, offset, rad1, rad2)
     Polygon *P;
     int n;
     float ctrx, ctry, offset, rad1, rad2;
{
  register int i;
  
  P->n = n;
  for (i = 0; i < n; i++)
    {
      Pv(i)[XX] = ctrx + rad1 * cos(offset + (TWO_PI * i) / n);
      Pv(i)[YY] = ctry + rad2 * sin(offset + (TWO_PI * i) / n);
    }
}


void Parse(argc, argv)
     int argc;
     char **argv;
{
  --argc; argv++;
  
  np = nq = 10;
  offp = offq = 0.0;
  rad1 = rad2 = rad3 =rad4 = 1.0;
  verbose = gr_verbose = randomize = 0;
  
  while ((argc > 0) && (argv[0][0] == '-'))
    {
      if (!strcmp(argv[0], "-n")) {
	if (sscanf(argv[1], "%d", &np) != 1) 
	  fatalError("ERROR -- \"%s\": bad n\n", argv[1]);
	if ((argc < 3) || (sscanf(argv[2], "%d", &nq) != 1))
	  {
	    nq = np;
	    argv += 2;
	    argc -= 2;
	  }
	else
	  {
	    argv += 3;
	    argc -= 3;
	  }
	if ((np > MAXPTS) || (nq > MAXPTS))
	  fatalError("ERROR -- \"%d\": bad n\n", (char *) nq);
	
      }
      else if (!strcmp(argv[0], "-o")) {
	if (sscanf(argv[1], "%f", &offp) != 1) 
	  fatalError("ERROR -- \"%s\": bad offset\n", argv[1]);
	if ((argc < 3) || (sscanf(argv[2], "%f", &offq) != 1)) 
	  {
	    offq = offp;
	    argv += 2;
	    argc -= 2;
	  }
	else
	  {
	    argv += 3;
	    argc -= 3;
	  }
      }
      else if (!strcmp(argv[0], "-R")) {
	if (sscanf(argv[1], "%f", &rad1) != 1)
	  fatalError("ERROR -- \"%s\": bad radii\n", argv[1]);
	if (sscanf(argv[2], "%f", &rad2) != 1)
	  fatalError("ERROR -- \"%s\": bad radii\n", argv[2]);
	if (sscanf(argv[3], "%f", &rad3) != 1)
	  fatalError("ERROR -- \"%s\": bad radii\n", argv[3]);
	if (sscanf(argv[4], "%f", &rad4) != 1)
	  fatalError("ERROR -- \"%s\": bad radii\n", argv[4]);
	argv += 5;
	argc -= 5;
      }
      else if (!strcmp(argv[0], "-g")) {
	if (verbose == 0) verbose = 1;
	if ((argc < 2) || (sscanf(argv[1], "%d", &gr_verbose) != 1))
	  gr_verbose = 1;
	else {
	  argv++;
	  --argc;
	}
	argv++;
	--argc;
      }
      else if (!strcmp(argv[0], "-v")) {
	if ((argc < 2) || (sscanf(argv[1], "%d", &verbose) != 1))
	  verbose++;
	else {
	  argv++;
	  --argc;
	}
	argv++;
	--argc;
      }
      else if (!strcmp(argv[0], "-r")) {
	argv++;
	--argc;
	randomize = 1;
      }
      else if (!strcmp(argv[0], "-h")) {
	argv++;
	--argc;
	fprintf(stderr, 
		"usage: main -n #P #Q -o circ_offsetP offQ -g{rafix} -v{erbose} -r{nd}\n"
		);
	exit(1);
      }
      
    }
}



void main(argc, argv)
     int argc;
     char **argv;
{
  HOMOG tangent;
  
  rad1 =  rad2 = rad3 = rad4 = 1.0;
  Parse(argc, argv);
  if (gr_verbose) 
    drawInit();
  do 
    {
      if (randomize)
	{
	  np = (random() & 0xFF) % (MAXPTS - 3) + 3 ;
	  nq = (random() & 0xFF) % (MAXPTS - 3) + 3;
	  offp = random() / 1.0e7;
	  offq = random() / 1.0e7;
	  rad1 = rsize[random() & 0x7];
	  rad2 = rsize[random() & 0x7];
	  rad3 = rsize[random() & 0x7];
	  rad4 = rsize[random() & 0x7];
	}
      getPoints(&P, np, 2.1, -0.05, offp, rad1, rad2);
      getPoints(&Q, nq, 0.0, 0.05, offq, rad3, rad4);
      
      printf("\nmain -n %4d %4d -o %g %g  -R %g %g %g %g", 
	     np, nq, offp, offq, rad1, rad2, rad3, rad4);
      fflush(stdout);
      Tang(&P, &Q);
      if (gr_verbose)
	{
	  color(WHITE);
	  drawSeg(Pvr(P.end), Qvr(Q.end));
	}
      displayAll();
      CROSSPROD_2CCH(Pvr(P.end), Qvr(Q.end), tangent); /* Get tangent */  
      if (!LEFT_PL(Pvr(P.end+1), tangent)) 
	{printf("\nPend+1 is wrong %g", 
		DOTPROD_2CH(Pvr(P.end+1), tangent)); randomize =0;} 
      if (RIGHT_PL(Pvr(P.end+P.n-1), tangent)) 
	{printf("\nPend-1 is wrong %g", 
		DOTPROD_2CH(Pvr(P.end+P.n-1), tangent)); randomize =0;}
      if (RIGHT_PL(Qvr(Q.end+1), tangent)) 
	{printf("\nQend+1 is wrong %g", 
		DOTPROD_2CH(Qvr(Q.end+1), tangent)); randomize =0;}
      if (!LEFT_PL(Qvr(Q.end+Q.n-1), tangent)) 
	{printf("\nQend-1 is wrong %g", 
		DOTPROD_2CH(Qvr(Q.end+Q.n-1), tangent)); randomize =0;}
      fflush(stdout);
      if (gr_verbose>1) sginap(25);/**/
    } while (randomize);
  sginap(500);
}

Makefile

# File Names
#CFLAGS = -O  -s 
CFLAGS = -g
LDFLAGS = -lc_s -lm -lgl_s
HDRS = main.h 
SRCS = main.c nosep.c
OBJS = main.o nosep.o
EXEC = main

all: $(OBJS) 
		cc $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(EXEC)

$(OBJS): $(HDRS)

.SUFFIXES:	.c .o .ln .h .out

.c.o:
		cc $(CFLAGS) -c $<