# 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

Include file defining data structure and mathematical operations

/* 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

The algorithm for computing a common tangent to two disjoint polygons

/* 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 */
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);

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

Calling program for debugging, display, & data generation.

/* 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 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 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);
}
}

#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;
{
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;
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)
if (sscanf(argv[2], "%f", &rad2) != 1)
if (sscanf(argv[3], "%f", &rad3) != 1)
if (sscanf(argv[4], "%f", &rad4) != 1)
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;

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",
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 \$<


-- Jack Snoeyink (snoeyink@cs.ubc.ca)