/* Generated by Interface Builder */

#include <stdio.h>
#include <stdlib.h> /* qsort blah blah */

#import "TwoDView.h"
#import "TwoDTSP.h"
#import <math.h>
#import <appkit/appkit.h>
#import <appkit/Matrix.h>
#import <appkit/Form.h>
#import <appkit/nextstd.h>
#import <dpsclient/wraps.h>

#define NI PSnewinstance()

#ifdef DEBUG /* put this in so the console sould not get flooded. */
// #define PRINTLOTS
#endif
#define MESSAGE(x) printf("%s\n",x)
@implementation TwoDView

/* Map the square [-1 .. 1]x[-1 .. 1] to the TwoDView, leaving a margin
 * of 10%, and centering the origin.
 *
 * the View must already be lockFocused.
 */
void setUpScaling(TwoDView *self)
{
  NXRect r;
  [self getBounds:&r];
  PStranslate(r.size.width/2,r.size.height/2);
  PSscale(0.9*0.5*r.size.width,0.9*0.5*r.size.height);
}

/*
 * set the drawing attributes (except for instance mode flag).
 * 
 * The code for setting line width averages the two sides of the
 * bounds rectangle, to avoid really nasty-looking lines when the view's
 * aspect ratio is not close to unity.  There are probably better ways
 * to do this.
 *
 * the View must already be lockFocused.
 */
void setUpDrawing(TwoDView *self)
{
  NXRect r;
  [self getBounds:&r];
  self->width=([self->lineWidthSlider floatValue]*2.0)/
              (r.size.width+r.size.height);
  if (self->width)
    PSsetlinewidth(self->width);
  else {
    if (NXDrawingStatus == NX_PRINTING) PSsetlinewidth(0.001);
    else PSsetlinewidth(0.0);
    };
}

/*
 * initialize instance variables relating to the data structures.
 * This is broken out of
 * the newFrame: method because mouseDown: events and open: will
 * also need to initialize those quantities.
 *
 * the instance variable numPoints contains the old value that
 * we will be replacing.
 *
 * cases: numPoints=npoints : no alloc, but clear flags.
 *        numPoints=0, npoints > 0 : alloc arrays, etc.
 *        numPoints>0, npoints = 0 : clobber arrays, etc.
 *        numPoints>0, npoints > 0 : realloc arrays, etc.
 *     
 * this function does NOT put any data in the array.  The caller must
 * do that.
 *
 * storage is not allocated for the Delaunay, Gabriel and relative neighborhood
 * graphs because they could be completely connected (large!) but usually
 * aren't.  They will be allocated in the appropriate doFooGraph routines
 * (below).
 */
void initGeometryStuff(TwoDView *self,int npoints)
{
#ifdef DEBUG
    MESSAGE("initGeometryStuff");
#endif

  if (self->numPoints == npoints) {
    self->mstExists=self->chExists=self->vorExists=
                    self->ggExists=self->rngExists=self->dtExists=NO;
    return;
    };
  if (!self->numPoints) { /* alloc */
    NX_MALLOC(self->data,float,npoints*2);
    NX_MALLOC(self->mst_edge,int,(npoints-1)*2);
    NX_MALLOC(self->ch_vertex,int,npoints);
    NX_MALLOC(self->tsp_edges,int,npoints*2);
    }
  else {
    if (!npoints) { /* clobber */
      NX_FREE(self->data);
      NX_FREE(self->mst_edge);
      NX_FREE(self->ch_vertex);
      NX_FREE(self->tsp_edges);
      }
    else {
      NX_REALLOC(self->data,float,npoints*2);
      NX_REALLOC(self->mst_edge,int,(npoints-1)*2);
      NX_REALLOC(self->tsp_edges,int,npoints*2);
      };
    };
  self->numPoints=npoints;
  [self->numPointsForm setIntValue:self->numPoints at:0];
  self->mstExists=self->chExists=self->vorExists=
                  self->ggExists=self->rngExists=self->dtExists=NO;
}

/*
 * make some text appear in the `Status' box
 */
#define STATUS(s) [statusText setStringValue:s]

/*
 * set things up.
 */
+ newFrame:(NXRect *)frameRect
{
  self = [[super newFrame:frameRect] allocateGState]; /* lots of focus pocus */
  numPoints=0; /* this will be fixed */
  operation=OP_PRIM_MST; /* agree with .nib */
  autoUpdate=1; /* agree with .nib */
  srandom(time(0));
  [self random:0];
  displayFlag = DISP_PTS_RESULTS;
  return self;
}

/* create a list of 2D points to demonstrate the algms on.
 * if sender is nil (zero), random: was called from initFrame. Otherwise,
 * the Generate button was pressed.
 */

- random:sender
{
  int i,n;

  if (sender && numPointsForm)
    n = [numPointsForm intValueAt:0];
  else
    n = 10; /* must be initializing */

  initGeometryStuff(self,n);  /* numPOints gets set */

  for(i=0;i<numPoints;i++) {
    X(i)=2.0*(random()&65535)/65536.0 - 1.0;
    Y(i)=2.0*(random()&65535)/65536.0 - 1.0;
    };

  if (sender) {
    if (autoUpdate) [self animThenDisp:0];
    else [self disp:0];
    };
  return self;
}

- mouseDown:(NXEvent *)theEvent
{
  NXPoint where = theEvent->location;
  NXRect r;
  float x,y,cx,cy;
  [self convertPoint:&where fromView:nil];
  [self getBounds:&r];
  cx=r.origin.x+0.5*r.size.width;
  cy=r.origin.y+0.5*r.size.height;

  x=(where.x - cx)/(0.5*0.9*r.size.width)+0.001;
  y=(where.y - cy)/(0.5*0.9*r.size.height)+0.001;
  if ((fabs(x) > 1.0)||(fabs(y) > 1.0)) return self;
#ifdef DEBUG
  fprintf(stderr,"x %f y %f\n",x,y);
#endif

  initGeometryStuff(self,numPoints+1); /* numPoints gets reset */
  X(numPoints-1)=x;
  Y(numPoints-1)=y;

  if (autoUpdate) [self animThenDisp:0];
  else [self disp:0];
  return self;
}

- free
{ initGeometryStuff(self,0); [super free]; return self; }
/*--------------------------------------------------------------------
Enable the form cell if it han an optimal value.
Assume that the cell is disabled.
                                                 Bill Roth/1991-Nov-26
--------------------------------------------------------------------*/

- setOperation:(Matrix *)sender
{ 
    operation=[sender selectedRow]; 
    if (operation >= OP_TSP_STUPID_NEIGHBOR) {
	[optimalValueFormCell setEnabled: YES];
	[timeSlider setEnabled: YES];
	[timeField setEnabled: YES];
    }
    else {
	[optimalValueFormCell setEnabled: NO];
	[timeSlider setEnabled: NO];
	[timeField setEnabled: NO];
	[optimizeButton setEnabled: NO];
    }
    return self; 
}

- setAutoUpdate:(Button *)sender
{ autoUpdate = [sender state]; return self; }


/* This method is invoked when the user bangs on the Animate/Display button. */
/*--------------------------------------------------------------------
Sets the display-pts-only flad, the calls [display] to draw points.
Then does the appropriate drawing routine, then when finished, calls 
[disp] which will call display and tell it do display everything.

                                                 Bill Roth/1991-Nov-17
--------------------------------------------------------------------*/

- animThenDisp:sender
{
#ifdef DEBUG
    MESSAGE("animthenDisp");
#endif

  /* erase prior drawing and display points only */
  displayFlag = DISP_PTS_ONLY;
  [self display];

  /* get ready to draw (in instance mode) */
  [self lockFocus];
  setUpScaling(self);
  NXSetColor([self->highColorWell color]);
  setUpDrawing(self);
  PSsetinstance(YES);

  /* do the operation (drawing in the doFoo methods is all instance drawing) */
  switch(operation) {
    case OP_PRIM_MST: [self doPrimMST]; break;
    case OP_KRUSKAL_MST: [self doKruskalMST]; break;
    case OP_JARVIS_HULL: [self doJarvisHull]; break;
    case OP_GRAHAM_HULL: [self doGrahamHull]; break;
    case OP_VORONOI_DIAGRAM: [self doVoronoiDiagram]; break;
    case OP_DELAUNAY_TRIANGULATION: [self doDelaunayTriangulation]; break;
    case OP_GABRIEL_GRAPH: [self doGabrielGraph]; break;
    case OP_RELATIVE_NEIGHBORHOOD_GRAPH: [self doRelativeNeighborhoodGraph];
                                         break;
    case OP_TSP_STUPID_NEIGHBOR:
			[self doSimpleNearNeighbor]; break;
    case OP_TSP_CHEAPEST_INSERTION: [self doCheapestInsertion]; break;
    case OP_TSP_NEAREST_NEIGHBOR: [self doNearestNeighbor]; break;
    case OP_TSP_NEAREST_ADDITION: [self doNearestAddition]; break;
    case OP_TSP_FARTHEST_INSERTION: [self doFarthestInsertion]; break;
    }

  /* clear instance mode drawing  and reset things*/
  PSsetinstance(NO);
  [self unlockFocus];

  /* display results */
  displayFlag = DISP_PTS_RESULTS;
  [self disp:0];
  return self;
}

-disp:sender
{
  [self display];
  return self;
}


/* Kruskal's MST algorithm (non-optimal quadratic version) */

static float *length;

int ycompar(p1,p2)
  int *p1,*p2;
  {
    float l1=length[*p1],l2=length[*p2];
    if (l1<l2) return -1;
    if (l1>l2) return 1;
    return 0;
  }

float cost(float *p1,float *p2) /* squd Eucl distance between 2 2D points */
{
  float x1,x2,y1,y2;
  x1= p1[0]; y1=p1[1];
  x2= p2[0]; y2=p2[1];
  return((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

/*
  XXX
       THIS IMPLEMENTATION IS NOT OPTIMAL
  XXX
*/

- doKruskalMST
{
  int *group,groupindex=1;
  int *edge,*edgptr,*index,idx;
  int i;
  float *lptr;
#ifdef DEBUG
    MESSAGE("doKruskalMST");
#endif

  STATUS("sorting edges by length");
  NX_MALLOC(edge,int,numPoints*(numPoints-1));
  edgptr=edge;
  NX_MALLOC(length,float,numPoints*(numPoints-1)/2);
  lptr=length;
  /* fill in edge endpoint and length arrays */
  for(i=0;i<numPoints-1;i++) {
    int j;
    char buf[80];
    sprintf(buf,"point %d",i);
    STATUS(buf);
    for(j=i+1;j<numPoints;j++) {
      *edgptr++ = i;
      *edgptr++ = j;
      *lptr++ = cost(XY(i),XY(j));
      };
    };
  STATUS("qsorting");
  NX_MALLOC(index,int,numPoints*(numPoints-1)/2);
  for(i=0;i<numPoints*(numPoints-1)/2;i++) index[i]=i;
  /* sort edges */
  qsort(index,numPoints*(numPoints-1)/2,sizeof(int),ycompar);

  NX_MALLOC(group,int,numPoints);
  bzero(group,numPoints*sizeof(int));
  
  /* Algorithm:
     FOR i = 1 to n-1 
       Add (to the tree) the shortest edge that does not form a circuit in
       the tree.
     ROF
  */
  idx=0;
  for(i=0;i<numPoints-1;i++) {
    char flag=FALSE;
    char buf[80];
    sprintf(buf,"adding edge %d",i);
    STATUS(buf);
    while (!flag) { /* look for an edge */
      int j;
      int v1=edge[index[idx]*2],g1=group[v1];
      int v2=edge[index[idx++]*2+1],g2=group[v2];
      NI;
      [self drawEdge:X(v1):Y(v1):X(v2):Y(v2)];
      NXPing();
      if ((!g1&&!g2)||(g1!=g2)) {
        flag=TRUE;
        mst_edge[i*2]=v1;
        mst_edge[i*2+1]=v2;
        if (g1&&!g2) group[v2]=g1;
        else if(g2&&!g1) group[v1]=g2;
        else if (!g1 && !g2) group[v1]=group[v2]=groupindex++;
        else  /* relabel all g2's to g1 */
          for(j=0;j<numPoints;j++) if (group[j]==g2) group[j]=g1;
        };
      };
    };
  NX_FREE(group);
  NX_FREE(index);
  NX_FREE(length);
  NX_FREE(edge);
  mstExists=YES;
  STATUS("idle");
  return self;
}


#define SIGNEDAREA(x1,y1,x2,y2,x3,y3) (x1*y2+x3*y1+x2*y3-x3*y2-x1*y3-x2*y1)

- doJarvisHull
{
/* Jarvis' march for planar convex hull - From Preparata and Shamos */
  float ymax= -3.0e30,ymin=3.0e30,y,best_x,best_y,candidate_x,candidate_y;
  float current_x,current_y,a;
  int i;
  int nh;
  int current_index=0,best_index=0;
  int min_index=0,max_index=0;
  char *used;
  /* 0. set up storage for indices */
#ifdef DEBUG
    MESSAGE("doJarvisHull");
#endif

  NX_MALLOC(used,char,numPoints);
  /* 1. find min & max in y-coordinate */
  STATUS("finding Ymin vertex");
  for(i=0;i<numPoints;i++) {
    y= Y(i);
    if (y<ymin) { ymin=y; min_index=i; };
    if (y>ymax) { ymax=y; max_index=i; };
    };
  nh = 1;

  ch_vertex[0]=min_index; /* first hull vertex */

  current_index=min_index;
  for(i=0;i<numPoints;i++) used[i]=FALSE;
  /* 2. Jarvis' march - min to max */ 
  STATUS("Jarvis' march: min to max");
  while (current_index!=max_index) {
    used[current_index]=TRUE;
    current_x= X(current_index);
    current_y= Y(current_index);
    best_x= current_x;
    best_y= current_y;
    for(i=0;i<numPoints;i++) if (!used[i]) {
      candidate_x= X(i);
      candidate_y= Y(i);
      NI;
      [self drawEdge:X(current_index):Y(current_index):X(i):Y(i)];
      NXPing();
      a=SIGNEDAREA(current_x,current_y,candidate_x,candidate_y, best_x,best_y);
      if (a>=0.0)
        { best_index=i; best_x=candidate_x; best_y=candidate_y;};
      };
    ch_vertex[nh++] = current_index = best_index;
    };
  /* 3. Jarvis' march - max to min */
  used[min_index]=FALSE;
  STATUS("Jarvis' march: max to min");
  while (current_index!=min_index) {
    used[current_index]=TRUE;
    current_x= X(current_index);
    current_y= Y(current_index);
    best_x=current_x;
    best_y=current_y;
    for(i=0;i<numPoints;i++)  if (!used[i]) {
      candidate_x= X(i);
      candidate_y= Y(i);
      NI;
      [self drawEdge:X(current_index):Y(current_index):candidate_x:candidate_y];
      NXPing();
      a=SIGNEDAREA(current_x,current_y,candidate_x,candidate_y, best_x,best_y);
      if (a>=0.0)
        { best_index=i; best_x=candidate_x; best_y=candidate_y;};
      };
    ch_vertex[nh++] = current_index = best_index;
    };
  nch_vertices=nh;
  NX_FREE(used);
  chExists=YES;
  STATUS("idle");
  return self;
}

/* Graham-hull utils */

#define QUAD(x,y) ((x<0) ? ((y<0)?3:2) : ((y<0)?4:1))

/* the following two declarations are nasty.  They are also the
 * easient way to solve a problem using qsort().  In the Jarvis
 * routine, the vertices are sorted by polar angle wrto the centroid
 * as origin.  However, I want to sort an array of indices to the
 * data rather than the data itself (for efficiency, etc.)  Since I need
 * access to the data by name, and the data is an instance variable, I
 * need to pass self to the plain-C routines and get the actual data by
 * following self.
 */
static float mx,my; /* used in xcompar and doJarvisHull */
static TwoDView *myself;

#define XX(i) myself->data[i*2]
#define YY(i) myself->data[i*2+1]

int xcompar(int *pp1,int *pp2) /*used in Graham scan to sort pts by polar angle*/
{
  int p1,p2,q1,q2;
  float x1,x2,y1,y2,a;
  p1= *pp1; p2= *pp2;
  x1= XX(p1)-mx; y1= YY(p1)-my;
  x2= XX(p2)-mx; y2= YY(p2)-my;
  /* cute li'l optimization.  If we can sort by quadrant, do it. */
  q1=QUAD(x1,y1); q2=QUAD(x2,y2);
  if (q1<q2)  return(-1);
  if (q2<q1)  return(1);
  /* if we get here, p1 & p2 are in same quadrant. shucks. */
  a=SIGNEDAREA(0.0,0.0,x2,y2,x1,y1);
  if (a<0.0) return(-1);
  if (a==0.0) return(0);
  return(1);
}

/* be safe. */
#undef XX
#undef YY

- doGrahamHull
{
  /* Graham's scan for planar convex hull. From Preparata and Shamos */
  int i,mindex=0,start,prev,current,next,nnext;
  int t1,t2,t3;
  float a,minx,xx;
  int *pointers;
  int nh;

  NXSetColor([highColorWell color]);

  /* THIS IS GROSS. */
  myself = self;

  NX_MALLOC(pointers,int,numPoints);
  nh=numPoints;
  /* find point in hull interior (centroid will do)*/
  STATUS("finding interior point");
  mx=0.0; my=0.0; minx=999;
  for(i=0;i<numPoints;i++) {
    mx += X(i);
    my += Y(i);
    };
  mx /= numPoints;
  my /= numPoints;

  for(i=0;i<numPoints;i++) {
    xx= X(i);
    if (xx<minx) { minx=xx; mindex=i;};
    };

  STATUS("sorting points by polar angle");
  /* sort points by polar angle and distance */
  for(i=0;i<numPoints;i++)  pointers[i]=i;
  qsort((void *)pointers,(size_t)numPoints,sizeof(int),xcompar);

  /* find min-X vertex (guaranteed to be on CH) */
  for(i=0;i<numPoints;i++) if (pointers[i]==mindex) break;

  start=current=i; /* point i is a hull vertex. */

  /* set up pointers in linked list */
  prev=(current==0)?nh-1:current-1;
  next=(current==nh-1)?0:current+1;
  nnext=(next==nh-1)?0:next+1;
  /* loop over all points */
  STATUS("Performing Graham's scan");
  while(pointers[next]!=pointers[start]) {
    if (nh==3) break; /* triangle is ALWAYS a convex polygon */
    t1= pointers[current]; t2= pointers[next]; t3= pointers[nnext];
    NI;
    [self drawEdge:X(current):Y(current):X(next):Y(next)];
    [self drawEdge:X(current):Y(current):X(nnext):Y(nnext)];
    [self drawEdge:X(nnext):Y(nnext):X(next):Y(next)];
    NXPing();
    a=SIGNEDAREA(X(t1),Y(t1),X(t2),Y(t2),X(t3),Y(t3));
    if (a<0.0) { /* delete next vertex */
      for(i=next+1;i<nh;i++) pointers[i-1]= pointers[i];
      nh--;
      if (start>=next) start--;
      if (current>=next) current--;
      if (current != start) {
        /* back up one vertex in list */
        current --;
        if (current<0) current += nh;
        prev=(current==0)?nh-1:current-1;
        next=(current==nh-1)?0:current+1;
        nnext=(next==nh-1)?0:next+1;
        }
      else { /* do not back up if current vertex is starting vertex. */
        prev=(current==0)?nh-1:current-1;
        next=(current==nh-1)?0:current+1;
        nnext=(next==nh-1)?0:next+1;
        };
      }
    else { /* move forward in linked list */
      prev=current;
      current=next;
      next=(next==nh-1)?0:next+1;
      nnext=(next==nh-1)?0:next+1;
      };
    };
  for(i=0;i<nh;i++) ch_vertex[i]=pointers[i];
  nch_vertices=nh;
  NX_FREE(pointers);
  chExists=YES;
  return self;
}

- doPrimMST
{
  /* Prim's quadratic MST algorithm - from Horowitz & Sahni */
  int *near;
  int i,j=0,k=0,l=0;
  float mincost,x,xmin;

#ifdef DEBUG
    MESSAGE("doPrimMST");
#endif

  NX_MALLOC(near,int,numPoints);
  bzero(near,numPoints*sizeof(int));
  mincost=10000000000.0;
  /* find shortest edge (O(n^2 time)) */
  STATUS("finding shortest edge");
  for(i=0;i<numPoints-1;i++) {
   char buf[80];
   sprintf(buf,"point %d",i);
   STATUS(buf);
   for(j=i+1;j<numPoints;j++) {
     x=cost(XY(i),XY(j));
     NI;
     [self drawEdge:X(i):Y(i):X(j):Y(j)];
     NXPing();
     if (x<mincost) { mincost=x; k=i; l=j; };
     };
   };
  /* (k,l) is the first edge in the MST */
  mst_edge[0]=k; mst_edge[1]=l;
  NI;
  [self drawEdge:X(k):Y(k):X(l):Y(l)];
  /* set up initial nearest-neighbor array for k and l */
  for(i=0;i<numPoints;i++)
    if (cost(XY(i),XY(l))<cost(XY(i),XY(k)))
      near[i]=l;
    else
      near[i]=k;
  near[k]= near[l]= -1;
  /* add the remaining n-2 edges */
  for(i=1;i<numPoints-1;i++) {
    char buf[80];
    sprintf(buf,"adding edge %d",i);
    STATUS(buf);
    xmin=100000000000.0;
    for(k=0;k<numPoints;k++) 
      if (near[k]!= -1) {
        x=cost(XY(k),XY(near[k]));
        if (x<xmin) { xmin=x; j=k;};
        };
    mst_edge[i*2] = j; mst_edge[i*2+1]= near[j];
    NI;
    [self drawEdge:X(j):Y(j):X(near[j]):Y(near[j])];
    mincost += xmin;
    near[j]= -1;
    for(k=0;k<numPoints;k++)
      if (near[k]!= -1)
        if (cost(XY(k),XY(near[k])) > cost(XY(k),XY(j)))
          near[k]=j;
    };
  NX_FREE(near);
  NXPing();
  mstExists=YES;
  STATUS("idle");
  return self;
}

- erase
{
  NXRect r;
  [[self getBounds:&r] lockFocus];
  NXSetColor([bgColorWell color]);
  NXRectFill(&r);
  [self unlockFocus];
  return self;
}

- drawDT  /* focus must already be locked */
{
  int i;
  for(i=0;i<ndt_edge;i++) {
    int i1=dt_edge[i*2],i2=dt_edge[i*2+1];
    char buf[80];
    sprintf(buf,"DT edge %d",i);
    STATUS(buf);
    [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
    };
  return self;
}

- drawGG  /* focus must already be locked */
{
  int i;
  for(i=0;i<ngg_edge;i++) {
    int i1=gg_edge[i*2],i2=gg_edge[i*2+1];
    char buf[80];
    sprintf(buf,"GG edge %d",i);
    STATUS(buf);
    [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
    };
  return self;
}

- drawRNG  /* focus must already be locked */
{
  int i;
  for(i=0;i<nrng_edge;i++) {
    int i1=rng_edge[i*2],i2=rng_edge[i*2+1];
    char buf[80];
    sprintf(buf,"RNG edge %d",i);
    STATUS(buf);
    [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
    };
  return self;
}

/*--------------------------------------------------------------------
This function does the actual drawing for the algorithms. Most of the 
routines merely draw line between sets of vertices.

                                                 Bill Roth/1991-Nov-17
--------------------------------------------------------------------*/
- drawSelf:(NXRect *)rects:(int)rectCount
{
  int i;
#ifdef DEBUG
  fprintf(stderr,"drawSelf (%d points)\n",numPoints);
#endif

  STATUS("redisplaying");
  NXSetColor([bgColorWell color]);
  NXRectFill(&rects[0]); //fill page with current color.

  setUpScaling(self); // set up scaling

  setUpDrawing(self);
  NXSetColor([self->fgColorWell color]);
  for(i=0;i<numPoints;i++) {
#ifdef PRINTLOTS
    fprintf(stderr,"%d: %f %f\n",i,X(i),Y(i));
#endif
    PSnewpath();
    PSarc(X(i),Y(i),width,0.0,360.0);
    PSfill();
    };
  NXPing();

  /* draw selected data structure */
  if (displayFlag != DISP_PTS_ONLY) {
#ifdef DEBUG
    MESSAGE("drawing everything");
#endif

    switch(operation) {
      case OP_PRIM_MST:
      case OP_KRUSKAL_MST: if (mstExists) {
        for(i=0;i<numPoints-1;i++) {
          int i1=mst_edge[i*2], i2=mst_edge[i*2+1];
          PSnewpath();
          PSmoveto(X(i1),Y(i1));
          PSlineto(X(i2),Y(i2));
          PSstroke();
          NXPing();
          };
        break;
        };
      case OP_JARVIS_HULL:
      case OP_GRAHAM_HULL: if (chExists) {
        PSnewpath();
        PSmoveto(X(ch_vertex[0]),Y(ch_vertex[0]));
        for(i=1;i<nch_vertices;i++) PSlineto(X(ch_vertex[i]),Y(ch_vertex[i]));
        PSclosepath();
        PSstroke();
        NXPing();
        break;
        };
      case OP_VORONOI_DIAGRAM: if (vorExists) {
        float *foo;
        NX_MALLOC(foo,float,numPoints*2);
        for(i=0;i<numPoints;i++) {
          int nv,j;
          if ((nv=find_vregion(i,foo)) == -1) {
            NXRunAlertPanel("Voronoi","find_vregion (%d) failed",
                            NULL,NULL,NULL,i);
            break;
            };
          PSnewpath();
          PSmoveto(foo[0],foo[1]);
          for(j=1;j<nv;j++) PSlineto(foo[j*2],foo[j*2+1]);
          PSclosepath();
          PSstroke();
          NXPing();
          };
        NX_FREE(foo);
        break;
        };
      case OP_DELAUNAY_TRIANGULATION: if (vorExists) [self drawDT]; break;
      case OP_GABRIEL_GRAPH: if (ggExists)  [self drawGG]; break;
      case OP_RELATIVE_NEIGHBORHOOD_GRAPH: if (rngExists) [self drawRNG]; break;
      case OP_TSP_STUPID_NEIGHBOR:
      case OP_TSP_NEAREST_NEIGHBOR: 
      case OP_TSP_CHEAPEST_INSERTION:
      case OP_TSP_NEAREST_ADDITION:
      case OP_TSP_FARTHEST_INSERTION:
	 [self drawTSP]; break;
      };
    };
  NXPing();
  STATUS("idle");
  return self;
}

/* draw an edge with given endpoints */
-drawEdge:(float)x1:(float)y1:(float)x2:(float)y2
{ /*assumes focus is locked already and attributes (color, width) already set*/
  PSnewpath();
  PSmoveto(x1,y1);
  PSlineto(x2,y2);
  PSstroke();
  return self;
}

- doVoronoiDiagram
{
  NXRect r;
  float xmin,xmax,ymin,ymax;

  STATUS("Finding Voronoi diagram");
  [self getBounds:&r];
  xmin=-1.01;
  xmax=1.01;
  ymin=-1.01;
  ymax=1.01;

  if (load_vsites(numPoints,data,xmin,ymin,xmax,ymax)) {
    NXRunAlertPanel("Voronoi","load_vsites failed",NULL,NULL,NULL);
    vorExists=NO;
    return self;
    };
  vorExists=YES;
  return self;
}

static void checkedge(TwoDView *s,int i1,int i2)
{
  int i;
  int ii1=MIN(i1,i2),ii2=MAX(i1,i2);
  for(i=0;i<s->ndt_edge;i++)
    if ((s->dt_edge[i*2]==ii1)&&(s->dt_edge[i*2+1]==ii2)) return;
  s->dt_edge[s->ndt_edge*2]=ii1;
  s->dt_edge[s->ndt_edge*2+1]=ii2;
  s->ndt_edge++;
}

- doDelaunayTriangulation
{
  int ndt,i;
  TRI *tris;
  [self doVoronoiDiagram]; /* cheat! */

  STATUS("building Delaunay triangulation data structures");
  /* from list of triangles, obtain minimal list of edged (toss dups) */
  if ((ndt=find_dtriangles(&tris)) == -1) {
    NXRunAlertPanel("Delaunay","find_dtriangles() failed",
                    NULL,NULL,NULL);
    return self;
    };

  NX_MALLOC(dt_edge,int,2*3*ndt);
  ndt_edge=0;

  for(i=0;i<ndt;i++) {
    int p1=tris[i].v1->pid;
    int p2=tris[i].v2->pid;
    int p3=tris[i].v3->pid;
    checkedge(self,p1,p2);
    checkedge(self,p1,p3);
    checkedge(self,p2,p3);
    };
  dtExists=YES;
#ifdef PRINTLOTS
  for(i=0;i<ndt_edge;i++)
    fprintf(stderr,"%d: %d %d\n",i,dt_edge[i*2],dt_edge[i*2+1]);
#endif
  return self;
}


- doGabrielGraph
{
  int i;

  if (!dtExists)[self doDelaunayTriangulation];
  /* we need to see the Delaunay graph while we're animating */
  NXSetColor([fgColorWell color]);
  PSsetinstance(NO);
  [[[self drawDT] window] flushWindow];
  NXSetColor([highColorWell color]);
  PSsetinstance(YES);

  ngg_edge=ndt_edge;
  NX_MALLOC(gg_edge,int,2*ngg_edge);

  bcopy(dt_edge,gg_edge,2*ngg_edge*sizeof(int));


  /*for each DT edge, check circle of influence */
  for(i=0;i<ngg_edge;i++) {
    int i1=gg_edge[i*2], i2=gg_edge[i*2+1];
    int j;
    float x1=X(i1),y1=Y(i1);
    float x2=X(i2),y2=Y(i2);
    float cx=(x1+x2)/2,cy=(y1+y2)/2,r=0.5*sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
    char buf[80];
#ifdef PRINTLOTS
    fprintf(stderr,"edge %d (%d %d)\n",i,i1,i2);
    fprintf(stderr,"(%f %f) -> (%f %f)\n",x1,y2,x2,y2);
    fprintf(stderr,"center %f %f rad %f\n",cx,cy,r);
#endif
    sprintf(buf,"considering Delaunay edge %d",i);
    STATUS(buf);
    NI;
    PSnewpath();
    PSarc(cx,cy,r,0.0,360.0);
    PSstroke();
    NXPing();

    /* check to see if any other vertex in within the circle of influence.
       if so, remove the edge from the GG */
    for(j=0;j<numPoints;j++) {
      if ((j!=i1)&&(j!=i2)) {
        float x=X(j),y=Y(j),d=(x-cx)*(x-cx)+(y-cy)*(y-cy)-r*r;
#ifdef PRINTLOTS
        fprintf(stderr,"checking %d (%f %f) d %g\n",j,x,y,d);
#endif
        if (d<0.0) { /* clobber edge */
          int k;
#ifdef DEBUG
          fprintf(stderr,"die.\n");
#endif
          for(k=i;k<ngg_edge-1;k++) {
            gg_edge[k*2]=gg_edge[(k+1)*2];
            gg_edge[k*2+1]=gg_edge[(k+1)*2+1];
            };
          ngg_edge--;
          /* highlight the dead edge */
          PSsetinstance(NO);
          PSnewpath();
          PSmoveto(x1,y1);
          PSlineto(x2,y2);
          PSstroke();
          PSsetinstance(YES);
          i--;
          break;
          };
        };
      };
    };
  NX_REALLOC(gg_edge,int,2*ngg_edge);
  ggExists=YES;
  return self;
}

- doRelativeNeighborhoodGraph
{
  int i;

  if (!dtExists)[self doDelaunayTriangulation];
  /* we need to see the Delaunay graph while we're animating */
  PSsetinstance(NO);
  NXSetColor([fgColorWell color]);
  [[[self drawDT] window] flushWindow];
  NXSetColor([highColorWell color]);
  PSsetinstance(YES);

  nrng_edge=ndt_edge;
  NX_MALLOC(rng_edge,int,2*nrng_edge);

  bcopy(dt_edge,rng_edge,2*nrng_edge*sizeof(int));


  /*for each DT edge, check lune of influence */
  for(i=0;i<nrng_edge;i++) {
    int i1=rng_edge[i*2], i2=rng_edge[i*2+1];
    int j;
    float x1=X(i1),y1=Y(i1);
    float x2=X(i2),y2=Y(i2);
    float r2=(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);
    float r=sqrt(r2);
    char buf[80];
#ifdef PRINTLOTS
    fprintf(stderr,"edge %d (%d %d)\n",i,i1,i2);
    fprintf(stderr,"(%f %f) -> (%f %f)\n",x1,y2,x2,y2);
    fprintf(stderr,"rad^2 %f\n",r2);
#endif
    sprintf(buf,"considering Delaunay edge %d",i);
    STATUS(buf);
    NI;
    PSnewpath();
    PSarc(x1,y1,r,0.0,360.0);
    PSmoveto(x2,y2); /* eliminate annoying connector */
    PSarc(x2,y2,r,0.0,360.0);
    PSstroke();
    NXPing();

    /* check to see if any other vertex in within the lune of influence.
       if so, remove the edge from the RNG */
    for(j=0;j<numPoints;j++) {
      if ((j!=i1)&&(j!=i2)) {
        float x=X(j),y=Y(j);
        float d1=(x-x1)*(x-x1)+(y-y1)*(y-y1)-r*r;
        float d2=(x-x2)*(x-x2)+(y-y2)*(y-y2)-r*r;
#ifdef PRINTLOTS
        fprintf(stderr,"checking %d (%f %f) d1 %g d2 %g\n",j,x,y,d1,d2);
#endif
        if ((d1<0.0)&&(d2<0.0)) { /* clobber edge */
          int k;
#ifdef PRINTLOTS
          fprintf(stderr,"die.\n");
#endif
          for(k=i;k<nrng_edge-1;k++) {
            rng_edge[k*2]=rng_edge[(k+1)*2];
            rng_edge[k*2+1]=rng_edge[(k+1)*2+1];
            };
          nrng_edge--;
          /* highlight the dead edge */
          PSsetinstance(NO);
          PSnewpath();
          PSmoveto(x1,y1);
          PSlineto(x2,y2);
          PSstroke();
          PSsetinstance(YES);
          i--;
          break;
          };
        };
      };
    };
  NX_REALLOC(rng_edge,int,2*nrng_edge);
  rngExists=YES;
  return self;
}

- saveTIFF:sender
{
  NXRect r;
  id bm=[NXBitmapImageRep alloc];
  id sp=[SavePanel new];
  NXStream *stream;
  int fd;

  [self lockFocus];
  [self getBounds:&r];
  [bm initData:NULL fromRect:&r];
  [self unlockFocus];
  
  [sp setRequiredFileType:"tiff"];
  if (![sp runModal]) return self;
  fd=open([sp filename],O_CREAT|O_RDWR,0644);
  /*[sp free]; screaming death! screaming death! don't do this! */
  stream = NXOpenFile(fd,NX_READWRITE);
  if (!stream) {
    NXRunAlertPanel("TIFF","NXOpenFile failed",NULL,NULL,NULL);
    perror("NXOpenFile");
    [bm free];
    close(fd);
    return self;
    };
  STATUS("saving image as TIFF");
  [bm writeTIFF:stream usingCompression:NX_TIFF_COMPRESSION_LZW];
  [bm free];
  NXClose(stream);
  close(fd);
  STATUS("idle");
  return self;
}

- saveData:sender
{
  id sp=[SavePanel new];
  FILE *fp;
  int i;

  [sp setRequiredFileType:"2Ddata"];
  if (![sp runModal]) return self;
  fp=fopen([sp filename],"w");
  /*[sp free];*/
  if (!fp) {
    NXRunAlertPanel("Data","Couldn't open file",NULL,NULL,NULL);
    return self;
    };
  STATUS("saving data");
  fprintf(fp,"%d\n",numPoints);
  for(i=0;i<numPoints;i++) fprintf(fp,"%g %g\n",X(i),Y(i));
  fclose(fp);
  STATUS("idle");
  return self;
}

- openData:sender
{
  id op=[OpenPanel new];
  const char *const types[] = { "2Ddata" , NULL};
  FILE *fp;
  int i;

  [op setRequiredFileType:"2Ddata"];
  [op runModalForDirectory:"." file:"" types:types];
  chdir([op directory]);
  if (![op filenames]) {
    /*[op free];*/
    return self;
    };

  fp=fopen(*[op filenames],"r");
  /*[op free];*/
  if (!fp) {
    NXRunAlertPanel("Open","Couldn't open file",NULL,NULL,NULL);
    return self;
    };

  STATUS("reading data");
  fscanf(fp,"%d ",&numPoints);
  initGeometryStuff(self,numPoints);
  for(i=0;i<numPoints;i++) {
    fscanf(fp,"%f %f ",&X(i),&Y(i));
    };
  fclose(fp);
  [self disp:0];
  STATUS("idle");
  return self;
}

- close:sender
{
  initGeometryStuff(self,0);
  [self erase];
  return self;
}
  
-lineWidthChanged:sender
{
  char buf[80];
  sprintf(buf,"Line width: %6.3f\n",[sender floatValue]);
  [lineWidthText setStringValue:buf];
  /* could redraw self here */
  if (autoUpdate) [self display];
  return self;
}


-setFgColorWell:sender
{
  fgColorWell = sender;
  [sender setColor:NXConvertGrayToColor(0.0)];
  return self;
}

-setBgColorWell:sender
{
  bgColorWell = sender;
  [sender setColor:NXConvertGrayToColor(1.0)];
  return self;
}


-setHighColorWell:sender
{
  highColorWell = sender;
  [sender setColor:NXConvertGrayToColor(0.666)];
  return self;
}

- (BOOL)acceptsFirstResponder { return YES; }
- (BOOL)acceptsFirstMouse { return YES; }

- copy:sender
{
  id pb=[Pasteboard new];
  NXStream *stream;
  char *d;
  int length,maxLength;

  STATUS("copying image to clipboard");
  [pb declareTypes:&NXPostScriptPboard num:1 owner:self];
  stream=NXOpenMemory(NULL,0,NX_WRITEONLY);
  [self copyPSCodeInside:NULL to:stream];
  NXGetMemoryBuffer(stream,&d,&length,&maxLength);
  [pb writeType:NXPostScriptPboard data:d length:length];
  NXCloseMemory(stream,NX_FREEBUFFER);
  [pb free];
  STATUS("idle");
  return self;
}

- saveEPS:sender
{
  NXStream *stream;
  id sp = [SavePanel new];
  int fd;

  [sp setRequiredFileType:"eps"];
  if (![sp runModal]) return self;
  fprintf(stderr,"filename %s\n",[sp filename]);
  fd=open([sp filename],O_CREAT|O_WRONLY,0644);
  /*[sp free];*/
  stream = NXOpenFile(fd,NX_WRITEONLY);
  if (!stream) {
    NXRunAlertPanel("EPS","NXOpenFile failed",NULL,NULL,NULL);
    perror("NXOpenFile");
    return self;
    };
  STATUS("saving image as EPS");
  [self copyPSCodeInside:NULL to:stream];
  NXClose(stream);
  close(fd);
  STATUS("idle");
  return self;
}
/*--------------------------------------------------------------------
Message from App object.

                                                 Bill Roth/1991-Dec-01
--------------------------------------------------------------------*/

- appDidInit: sender;
{
    drawTime = [timeSlider intValue];
    [timeField setIntValue: drawTime at: 0];
    [timeField setEnabled: NO]; // agree with NIB
	[timeSlider setEnabled: NO];// agree with NIB
    return self ;
}
- doOptimize :sender
{
	[self Optimize];
	return self;
}
@end
