home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Graphics / 2D_3D / 2DLab / Source / TwoDView.m < prev    next >
Encoding:
Text File  |  1995-06-12  |  28.3 KB  |  1,142 lines
  1. /* Generated by Interface Builder */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h> /* qsort blah blah */
  5.  
  6. #import "TwoDView.h"
  7. #import <math.h>
  8. #import <appkit/appkit.h>
  9. #import <appkit/Matrix.h>
  10. #import <appkit/Form.h>
  11. #import <appkit/nextstd.h>
  12. #import <dpsclient/wraps.h>
  13.  
  14. #define NI PSnewinstance()
  15.  
  16. @implementation TwoDView
  17.  
  18. /* Map the square [-1 .. 1]x[-1 .. 1] to the TwoDView, leaving a margin
  19.  * of 10%, and centering the origin.
  20.  *
  21.  * the View must already be lockFocused.
  22.  */
  23. void setUpScaling(TwoDView *self)
  24. {
  25.   NXRect r;
  26.   [self getBounds:&r];
  27.   PStranslate(r.size.width/2,r.size.height/2);
  28.   PSscale(0.9*0.5*r.size.width,0.9*0.5*r.size.height);
  29. }
  30.  
  31. /*
  32.  * set the drawing attributes (except for instance mode flag).
  33.  * 
  34.  * The code for setting line width averages the two sides of the
  35.  * bounds rectangle, to avoid really nasty-looking lines when the view's
  36.  * aspect ratio is not close to unity.  There are probably better ways
  37.  * to do this.
  38.  *
  39.  * the View must already be lockFocused.
  40.  */
  41. void setUpDrawing(TwoDView *self)
  42. {
  43.   NXRect r;
  44.   [self getBounds:&r];
  45.   self->width=([self->lineWidthSlider floatValue]*2.0)/
  46.               (r.size.width+r.size.height);
  47.   if (self->width)
  48.     PSsetlinewidth(self->width);
  49.   else {
  50.     if (NXDrawingStatus == NX_PRINTING) PSsetlinewidth(0.001);
  51.     else PSsetlinewidth(0.0);
  52.     };
  53. }
  54.  
  55. /*
  56.  * initialize instance variables relating to the data structures.
  57.  * This is broken out of
  58.  * the newFrame: method because mouseDown: events and open: will
  59.  * also need to initialize those quantities.
  60.  *
  61.  * the instance variable numPoints contains the old value that
  62.  * we will be replacing.
  63.  *
  64.  * cases: numPoints=npoints : no alloc, but clear flags.
  65.  *        numPoints=0, npoints > 0 : alloc arrays, etc.
  66.  *        numPoints>0, npoints = 0 : clobber arrays, etc.
  67.  *        numPoints>0, npoints > 0 : realloc arrays, etc.
  68.  *     
  69.  * this function does NOT put any data in the array.  The caller must
  70.  * do that.
  71.  *
  72.  * storage is not allocated for the Delaunay, Gabriel and relative neighborhood
  73.  * graphs because they could be completely connected (large!) but usually
  74.  * aren't.  They will be allocated in the appropriate doFooGraph routines
  75.  * (below).
  76.  */
  77. void initGeometryStuff(TwoDView *self,int npoints)
  78. {
  79.   if (self->numPoints == npoints) {
  80.     self->mstExists=self->chExists=self->vorExists=
  81.                     self->ggExists=self->rngExists=self->dtExists=NO;
  82.     return;
  83.     };
  84.   if (!self->numPoints) { /* alloc */
  85.     NX_MALLOC(self->data,float,npoints*2);
  86.     NX_MALLOC(self->mst_edge,int,(npoints-1)*2);
  87.     NX_MALLOC(self->ch_vertex,int,npoints);
  88.     }
  89.   else {
  90.     if (!npoints) { /* clobber */
  91.       NX_FREE(self->data);
  92.       NX_FREE(self->mst_edge);
  93.       NX_FREE(self->ch_vertex);
  94.       }
  95.     else {
  96.       NX_REALLOC(self->data,float,npoints*2);
  97.       NX_REALLOC(self->mst_edge,int,(npoints-1)*2);
  98.       NX_REALLOC(self->ch_vertex,int,npoints);
  99.       };
  100.     };
  101.   self->numPoints=npoints;
  102.   [self->numPointsForm setIntValue:self->numPoints at:0];
  103.   self->mstExists=self->chExists=self->vorExists=
  104.                   self->ggExists=self->rngExists=self->dtExists=NO;
  105. }
  106.  
  107. /*
  108.  * make some text appear in the `Status' box
  109.  */
  110. #define STATUS(s) [statusText setStringValue:s]
  111.  
  112. /*
  113.  * set things up.
  114.  */
  115. + newFrame:(NXRect *)frameRect
  116. {
  117.   self = [[super newFrame:frameRect] allocateGState]; /* lots of focus pocus */
  118.   numPoints=0; /* this will be fixed */
  119.   operation=OP_PRIM_MST; /* agree with .nib */
  120.   autoUpdate=1; /* agree with .nib */
  121.   srandom(time(0));
  122.   [self random:0];
  123.   displayFlag = DISP_PTS_RESULTS;
  124.   return self;
  125. }
  126.  
  127. /* create a list of 2D points to demonstrate the algms on.
  128.  * if sender is nil (zero), random: was called from initFrame. Otherwise,
  129.  * the Generate button was pressed.
  130.  */
  131.  
  132. - random:sender
  133. {
  134.   int i,n;
  135.  
  136.   if (sender && numPointsForm)
  137.     n = [numPointsForm intValueAt:0];
  138.   else
  139.     n = 10; /* must be initializing */
  140.  
  141.   initGeometryStuff(self,n);  /* numPOints gets set */
  142.  
  143.   for(i=0;i<numPoints;i++) {
  144.     X(i)=2.0*(random()&65535)/65536.0 - 1.0;
  145.     Y(i)=2.0*(random()&65535)/65536.0 - 1.0;
  146.     };
  147.  
  148.   if (sender) {
  149.     if (autoUpdate) [self animThenDisp:0];
  150.     else [self disp:0];
  151.     };
  152.   return self;
  153. }
  154.  
  155. - mouseDown:(NXEvent *)theEvent
  156. {
  157.   NXPoint where = theEvent->location;
  158.   NXRect r;
  159.   float x,y,cx,cy;
  160.   [self convertPoint:&where fromView:nil];
  161.   [self getBounds:&r];
  162.   cx=r.origin.x+0.5*r.size.width;
  163.   cy=r.origin.y+0.5*r.size.height;
  164.  
  165.   x=(where.x - cx)/(0.5*0.9*r.size.width)+0.001;
  166.   y=(where.y - cy)/(0.5*0.9*r.size.height)+0.001;
  167.   if ((fabs(x) > 1.0)||(fabs(y) > 1.0)) return self;
  168. #ifdef DEBUG
  169.   fprintf(stderr,"x %f y %f\n",x,y);
  170. #endif
  171.  
  172.   initGeometryStuff(self,numPoints+1); /* numPoints gets reset */
  173.   X(numPoints-1)=x;
  174.   Y(numPoints-1)=y;
  175.  
  176.   if (autoUpdate) [self animThenDisp:0];
  177.   else [self disp:0];
  178.   return self;
  179. }
  180.  
  181. - free
  182. { initGeometryStuff(self,0); [super free]; return self; }
  183.  
  184. - setOperation:(Matrix *)sender
  185. { operation=[sender selectedRow]; return self; }
  186.  
  187. - setAutoUpdate:(Button *)sender
  188. { autoUpdate = [sender state]; return self; }
  189.  
  190.  
  191. /* This method is invoked when the user bangs on the Animate/Display button. */
  192.  
  193. - animThenDisp:sender
  194. {
  195.   /* erase prior drawing and display points only */
  196.   displayFlag = DISP_PTS_ONLY;
  197.   [self display];
  198.  
  199.   /* get ready to draw (in instance mode) */
  200.   [self lockFocus];
  201.   setUpScaling(self);
  202.   NXSetColor([self->highColorWell color]);
  203.   setUpDrawing(self);
  204.   PSsetinstance(YES);
  205.  
  206.   /* do the operation (drawing in the doFoo methods is all instance drawing) */
  207.   switch(operation) {
  208.     case OP_PRIM_MST: [self doPrimMST]; break;
  209.     case OP_KRUSKAL_MST: [self doKruskalMST]; break;
  210.     case OP_JARVIS_HULL: [self doJarvisHull]; break;
  211.     case OP_GRAHAM_HULL: [self doGrahamHull]; break;
  212.     case OP_VORONOI_DIAGRAM: [self doVoronoiDiagram]; break;
  213.     case OP_DELAUNAY_TRIANGULATION: [self doDelaunayTriangulation]; break;
  214.     case OP_GABRIEL_GRAPH: [self doGabrielGraph]; break;
  215.     case OP_RELATIVE_NEIGHBORHOOD_GRAPH: [self doRelativeNeighborhoodGraph];
  216.                                          break;
  217.     }
  218.  
  219.   /* clear instance mode drawing  and reset things*/
  220.   PSsetinstance(NO);
  221.   [self unlockFocus];
  222.  
  223.   /* display results */
  224.   displayFlag = DISP_PTS_RESULTS;
  225.   [self disp:0];
  226.   return self;
  227. }
  228.  
  229. -disp:sender
  230. {
  231.   [self display];
  232.   return self;
  233. }
  234.  
  235.  
  236. /* Kruskal's MST algorithm (non-optimal quadratic version) */
  237.  
  238. static float *length;
  239.  
  240. int ycompar(p1,p2)
  241.   int *p1,*p2;
  242.   {
  243.     float l1=length[*p1],l2=length[*p2];
  244.     if (l1<l2) return -1;
  245.     if (l1>l2) return 1;
  246.     return 0;
  247.   }
  248.  
  249. float cost(float *p1,float *p2) /* squd Eucl distance between 2 2D points */
  250. {
  251.   float x1,x2,y1,y2;
  252.   x1= p1[0]; y1=p1[1];
  253.   x2= p2[0]; y2=p2[1];
  254.   return((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  255. }
  256.  
  257. /*
  258.   XXX
  259.        THIS IMPLEMENTATION IS NOT OPTIMAL
  260.   XXX
  261. */
  262.  
  263. - doKruskalMST
  264. {
  265.   int *group,groupindex=1;
  266.   int *edge,*edgptr,*index,idx;
  267.   int i;
  268.   float *lptr;
  269.  
  270.   STATUS("sorting edges by length");
  271.   NX_MALLOC(edge,int,numPoints*(numPoints-1));
  272.   edgptr=edge;
  273.   NX_MALLOC(length,float,numPoints*(numPoints-1)/2);
  274.   lptr=length;
  275.   /* fill in edge endpoint and length arrays */
  276.   for(i=0;i<numPoints-1;i++) {
  277.     int j;
  278.     char buf[80];
  279.     sprintf(buf,"point %d",i);
  280.     STATUS(buf);
  281.     for(j=i+1;j<numPoints;j++) {
  282.       *edgptr++ = i;
  283.       *edgptr++ = j;
  284.       *lptr++ = cost(XY(i),XY(j));
  285.       };
  286.     };
  287.   STATUS("qsorting");
  288.   NX_MALLOC(index,int,numPoints*(numPoints-1)/2);
  289.   for(i=0;i<numPoints*(numPoints-1)/2;i++) index[i]=i;
  290.   /* sort edges */
  291.   qsort(index,numPoints*(numPoints-1)/2,sizeof(int),ycompar);
  292.  
  293.   NX_MALLOC(group,int,numPoints);
  294.   bzero(group,numPoints*sizeof(int));
  295.   
  296.   /* Algorithm:
  297.      FOR i = 1 to n-1 
  298.        Add (to the tree) the shortest edge that does not form a circuit in
  299.        the tree.
  300.      ROF
  301.   */
  302.   idx=0;
  303.   for(i=0;i<numPoints-1;i++) {
  304.     char flag=FALSE;
  305.     char buf[80];
  306.     sprintf(buf,"adding edge %d",i);
  307.     STATUS(buf);
  308.     while (!flag) { /* look for an edge */
  309.       int j;
  310.       int v1=edge[index[idx]*2],g1=group[v1];
  311.       int v2=edge[index[idx++]*2+1],g2=group[v2];
  312.       NI;
  313.       [self drawEdge:X(v1):Y(v1):X(v2):Y(v2)];
  314.       NXPing();
  315.       if ((!g1&&!g2)||(g1!=g2)) {
  316.         flag=TRUE;
  317.         mst_edge[i*2]=v1;
  318.         mst_edge[i*2+1]=v2;
  319.         if (g1&&!g2) group[v2]=g1;
  320.         else if(g2&&!g1) group[v1]=g2;
  321.         else if (!g1 && !g2) group[v1]=group[v2]=groupindex++;
  322.         else  /* relabel all g2's to g1 */
  323.           for(j=0;j<numPoints;j++) if (group[j]==g2) group[j]=g1;
  324.         };
  325.       };
  326.     };
  327.   NX_FREE(group);
  328.   NX_FREE(index);
  329.   NX_FREE(length);
  330.   NX_FREE(edge);
  331.   mstExists=YES;
  332.   STATUS("idle");
  333.   return self;
  334. }
  335.  
  336.  
  337. #define SIGNEDAREA(x1,y1,x2,y2,x3,y3) (x1*y2+x3*y1+x2*y3-x3*y2-x1*y3-x2*y1)
  338.  
  339. - doJarvisHull
  340. {
  341. /* Jarvis' march for planar convex hull - From Preparata and Shamos */
  342.   float ymax= -3.0e30,ymin=3.0e30,y,best_x,best_y,candidate_x,candidate_y;
  343.   float current_x,current_y,a;
  344.   int i;
  345.   int nh;
  346.   int current_index=0,best_index=0;
  347.   int min_index=0,max_index=0;
  348.   char *used;
  349.   /* 0. set up storage for indices */
  350.   NX_MALLOC(used,char,numPoints);
  351.   /* 1. find min & max in y-coordinate */
  352.   STATUS("finding Ymin vertex");
  353.   for(i=0;i<numPoints;i++) {
  354.     y= Y(i);
  355.     if (y<ymin) { ymin=y; min_index=i; };
  356.     if (y>ymax) { ymax=y; max_index=i; };
  357.     };
  358.   nh = 1;
  359.  
  360.   ch_vertex[0]=min_index; /* first hull vertex */
  361.  
  362.   current_index=min_index;
  363.   for(i=0;i<numPoints;i++) used[i]=FALSE;
  364.   /* 2. Jarvis' march - min to max */ 
  365.   STATUS("Jarvis' march: min to max");
  366.   while (current_index!=max_index) {
  367.     used[current_index]=TRUE;
  368.     current_x= X(current_index);
  369.     current_y= Y(current_index);
  370.     best_x= current_x;
  371.     best_y= current_y;
  372.     for(i=0;i<numPoints;i++) if (!used[i]) {
  373.       candidate_x= X(i);
  374.       candidate_y= Y(i);
  375.       NI;
  376.       [self drawEdge:X(current_index):Y(current_index):X(i):Y(i)];
  377.       NXPing();
  378.       a=SIGNEDAREA(current_x,current_y,candidate_x,candidate_y, best_x,best_y);
  379.       if (a>=0.0)
  380.         { best_index=i; best_x=candidate_x; best_y=candidate_y;};
  381.       };
  382.     ch_vertex[nh++] = current_index = best_index;
  383.     };
  384.   /* 3. Jarvis' march - max to min */
  385.   used[min_index]=FALSE;
  386.   STATUS("Jarvis' march: max to min");
  387.   while (current_index!=min_index) {
  388.     used[current_index]=TRUE;
  389.     current_x= X(current_index);
  390.     current_y= Y(current_index);
  391.     best_x=current_x;
  392.     best_y=current_y;
  393.     for(i=0;i<numPoints;i++)  if (!used[i]) {
  394.       candidate_x= X(i);
  395.       candidate_y= Y(i);
  396.       NI;
  397.       [self drawEdge:X(current_index):Y(current_index):candidate_x:candidate_y];
  398.       NXPing();
  399.       a=SIGNEDAREA(current_x,current_y,candidate_x,candidate_y, best_x,best_y);
  400.       if (a>=0.0)
  401.         { best_index=i; best_x=candidate_x; best_y=candidate_y;};
  402.       };
  403.     ch_vertex[nh++] = current_index = best_index;
  404.     };
  405.   nch_vertices=nh;
  406.   NX_FREE(used);
  407.   chExists=YES;
  408.   STATUS("idle");
  409.   return self;
  410. }
  411.  
  412. /* Graham-hull utils */
  413.  
  414. #define QUAD(x,y) ((x<0) ? ((y<0)?3:2) : ((y<0)?4:1))
  415.  
  416. /* the following two declarations are nasty.  They are also the
  417.  * easient way to solve a problem using qsort().  In the Jarvis
  418.  * routine, the vertices are sorted by polar angle wrto the centroid
  419.  * as origin.  However, I want to sort an array of indices to the
  420.  * data rather than the data itself (for efficiency, etc.)  Since I need
  421.  * access to the data by name, and the data is an instance variable, I
  422.  * need to pass self to the plain-C routines and get the actual data by
  423.  * following self.
  424.  */
  425. static float mx,my; /* used in xcompar and doJarvisHull */
  426. static TwoDView *myself;
  427.  
  428. #define XX(i) myself->data[i*2]
  429. #define YY(i) myself->data[i*2+1]
  430.  
  431. int xcompar(int *pp1,int *pp2) /*used in Graham scan to sort pts by polar angle*/
  432. {
  433.   int p1,p2,q1,q2;
  434.   float x1,x2,y1,y2,a;
  435.   p1= *pp1; p2= *pp2;
  436.   x1= XX(p1)-mx; y1= YY(p1)-my;
  437.   x2= XX(p2)-mx; y2= YY(p2)-my;
  438.   /* cute li'l optimization.  If we can sort by quadrant, do it. */
  439.   q1=QUAD(x1,y1); q2=QUAD(x2,y2);
  440.   if (q1<q2)  return(-1);
  441.   if (q2<q1)  return(1);
  442.   /* if we get here, p1 & p2 are in same quadrant. shucks. */
  443.   a=SIGNEDAREA(0.0,0.0,x2,y2,x1,y1);
  444.   if (a<0.0) return(-1);
  445.   if (a==0.0) return(0);
  446.   return(1);
  447. }
  448.  
  449. /* be safe. */
  450. #undef XX
  451. #undef YY
  452.  
  453. - doGrahamHull
  454. {
  455.   /* Graham's scan for planar convex hull. From Preparata and Shamos */
  456.   int i,mindex=0,start,prev,current,next,nnext;
  457.   int t1,t2,t3;
  458.   float a,minx,xx;
  459.   int *pointers;
  460.   int nh;
  461.  
  462.   NXSetColor([highColorWell color]);
  463.  
  464.   /* THIS IS GROSS. */
  465.   myself = self;
  466.  
  467.   NX_MALLOC(pointers,int,numPoints);
  468.   nh=numPoints;
  469.   /* find point in hull interior (centroid will do)*/
  470.   STATUS("finding interior point");
  471.   mx=0.0; my=0.0; minx=999;
  472.   for(i=0;i<numPoints;i++) {
  473.     mx += X(i);
  474.     my += Y(i);
  475.     };
  476.   mx /= numPoints;
  477.   my /= numPoints;
  478.  
  479.   for(i=0;i<numPoints;i++) {
  480.     xx= X(i);
  481.     if (xx<minx) { minx=xx; mindex=i;};
  482.     };
  483.  
  484.   STATUS("sorting points by polar angle");
  485.   /* sort points by polar angle and distance */
  486.   for(i=0;i<numPoints;i++)  pointers[i]=i;
  487.   qsort((void *)pointers,(size_t)numPoints,sizeof(int),xcompar);
  488.  
  489.   /* find min-X vertex (guaranteed to be on CH) */
  490.   for(i=0;i<numPoints;i++) if (pointers[i]==mindex) break;
  491.  
  492.   start=current=i; /* point i is a hull vertex. */
  493.  
  494.   /* set up pointers in linked list */
  495.   prev=(current==0)?nh-1:current-1;
  496.   next=(current==nh-1)?0:current+1;
  497.   nnext=(next==nh-1)?0:next+1;
  498.   /* loop over all points */
  499.   STATUS("Performing Graham's scan");
  500.   while(pointers[next]!=pointers[start]) {
  501.     if (nh==3) break; /* triangle is ALWAYS a convex polygon */
  502.     t1= pointers[current]; t2= pointers[next]; t3= pointers[nnext];
  503.     NI;
  504.     [self drawEdge:X(current):Y(current):X(next):Y(next)];
  505.     [self drawEdge:X(current):Y(current):X(nnext):Y(nnext)];
  506.     [self drawEdge:X(nnext):Y(nnext):X(next):Y(next)];
  507.     NXPing();
  508.     a=SIGNEDAREA(X(t1),Y(t1),X(t2),Y(t2),X(t3),Y(t3));
  509.     if (a<0.0) { /* delete next vertex */
  510.       for(i=next+1;i<nh;i++) pointers[i-1]= pointers[i];
  511.       nh--;
  512.       if (start>=next) start--;
  513.       if (current>=next) current--;
  514.       if (current != start) {
  515.         /* back up one vertex in list */
  516.         current --;
  517.         if (current<0) current += nh;
  518.         prev=(current==0)?nh-1:current-1;
  519.         next=(current==nh-1)?0:current+1;
  520.         nnext=(next==nh-1)?0:next+1;
  521.         }
  522.       else { /* do not back up if current vertex is starting vertex. */
  523.         prev=(current==0)?nh-1:current-1;
  524.         next=(current==nh-1)?0:current+1;
  525.         nnext=(next==nh-1)?0:next+1;
  526.         };
  527.       }
  528.     else { /* move forward in linked list */
  529.       prev=current;
  530.       current=next;
  531.       next=(next==nh-1)?0:next+1;
  532.       nnext=(next==nh-1)?0:next+1;
  533.       };
  534.     };
  535.   for(i=0;i<nh;i++) ch_vertex[i]=pointers[i];
  536.   nch_vertices=nh;
  537.   NX_FREE(pointers);
  538.   chExists=YES;
  539.   return self;
  540. }
  541.  
  542. - doPrimMST
  543. {
  544.   /* Prim's quadratic MST algorithm - from Horowitz & Sahni */
  545.   int *near;
  546.   int i,j=0,k=0,l=0;
  547.   float mincost,x,xmin;
  548.   NX_MALLOC(near,int,numPoints);
  549.   bzero(near,numPoints*sizeof(int));
  550.   mincost=10000000000.0;
  551.   /* find shortest edge (O(n^2 time)) */
  552.   STATUS("finding shortest edge");
  553.   for(i=0;i<numPoints-1;i++) {
  554.    char buf[80];
  555.    sprintf(buf,"point %d",i);
  556.    STATUS(buf);
  557.    for(j=i+1;j<numPoints;j++) {
  558.      x=cost(XY(i),XY(j));
  559.      NI;
  560.      [self drawEdge:X(i):Y(i):X(j):Y(j)];
  561.      NXPing();
  562.      if (x<mincost) { mincost=x; k=i; l=j; };
  563.      };
  564.    };
  565.   /* (k,l) is the first edge in the MST */
  566.   mst_edge[0]=k; mst_edge[1]=l;
  567.   NI;
  568.   [self drawEdge:X(k):Y(k):X(l):Y(l)];
  569.   /* set up initial nearest-neighbor array for k and l */
  570.   for(i=0;i<numPoints;i++)
  571.     if (cost(XY(i),XY(l))<cost(XY(i),XY(k)))
  572.       near[i]=l;
  573.     else
  574.       near[i]=k;
  575.   near[k]= near[l]= -1;
  576.   /* add the remaining n-2 edges */
  577.   for(i=1;i<numPoints-1;i++) {
  578.     char buf[80];
  579.     sprintf(buf,"adding edge %d",i);
  580.     STATUS(buf);
  581.     xmin=100000000000.0;
  582.     for(k=0;k<numPoints;k++) 
  583.       if (near[k]!= -1) {
  584.         x=cost(XY(k),XY(near[k]));
  585.         if (x<xmin) { xmin=x; j=k;};
  586.         };
  587.     mst_edge[i*2] = j; mst_edge[i*2+1]= near[j];
  588.     NI;
  589.     [self drawEdge:X(j):Y(j):X(near[j]):Y(near[j])];
  590.     mincost += xmin;
  591.     near[j]= -1;
  592.     for(k=0;k<numPoints;k++)
  593.       if (near[k]!= -1)
  594.         if (cost(XY(k),XY(near[k])) > cost(XY(k),XY(j)))
  595.           near[k]=j;
  596.     };
  597.   NX_FREE(near);
  598.   NXPing();
  599.   mstExists=YES;
  600.   STATUS("idle");
  601.   return self;
  602. }
  603.  
  604. - erase
  605. {
  606.   NXRect r;
  607.   [[self getBounds:&r] lockFocus];
  608.   NXSetColor([bgColorWell color]);
  609.   NXRectFill(&r);
  610.   [self unlockFocus];
  611.   return self;
  612. }
  613.  
  614. - drawDT  /* focus must already be locked */
  615. {
  616.   int i;
  617.   for(i=0;i<ndt_edge;i++) {
  618.     int i1=dt_edge[i*2],i2=dt_edge[i*2+1];
  619.     char buf[80];
  620.     sprintf(buf,"DT edge %d",i);
  621.     STATUS(buf);
  622.     [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
  623.     };
  624.   return self;
  625. }
  626.  
  627. - drawGG  /* focus must already be locked */
  628. {
  629.   int i;
  630.   for(i=0;i<ngg_edge;i++) {
  631.     int i1=gg_edge[i*2],i2=gg_edge[i*2+1];
  632.     char buf[80];
  633.     sprintf(buf,"GG edge %d",i);
  634.     STATUS(buf);
  635.     [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
  636.     };
  637.   return self;
  638. }
  639.  
  640. - drawRNG  /* focus must already be locked */
  641. {
  642.   int i;
  643.   for(i=0;i<nrng_edge;i++) {
  644.     int i1=rng_edge[i*2],i2=rng_edge[i*2+1];
  645.     char buf[80];
  646.     sprintf(buf,"RNG edge %d",i);
  647.     STATUS(buf);
  648.     [self drawEdge:X(i1):Y(i1):X(i2):Y(i2)];
  649.     };
  650.   return self;
  651. }
  652.  
  653. - drawSelf:(NXRect *)rects:(int)rectCount
  654. {
  655.   int i;
  656. #ifdef DEBUG
  657.   fprintf(stderr,"drawSelf (%d points)\n",numPoints);
  658. #endif
  659.  
  660.   STATUS("redisplaying");
  661.   NXSetColor([bgColorWell color]);
  662.   NXRectFill(&rects[0]);
  663.  
  664.   setUpScaling(self);
  665.  
  666.   setUpDrawing(self);
  667.   NXSetColor([self->fgColorWell color]);
  668.   for(i=0;i<numPoints;i++) {
  669. #ifdef DEBUG
  670.     fprintf(stderr,"%d: %f %f\n",i,X(i),Y(i));
  671. #endif
  672.     PSnewpath();
  673.     PSarc(X(i),Y(i),width,0.0,360.0);
  674.     PSfill();
  675.     };
  676.   NXPing();
  677.  
  678.   /* draw selected data structure */
  679.   if (displayFlag != DISP_PTS_ONLY) {
  680.     switch(operation) {
  681.       case OP_PRIM_MST:
  682.       case OP_KRUSKAL_MST: if (mstExists) {
  683.         for(i=0;i<numPoints-1;i++) {
  684.           int i1=mst_edge[i*2], i2=mst_edge[i*2+1];
  685.           PSnewpath();
  686.           PSmoveto(X(i1),Y(i1));
  687.           PSlineto(X(i2),Y(i2));
  688.           PSstroke();
  689.           NXPing();
  690.           };
  691.         break;
  692.         };
  693.       case OP_JARVIS_HULL:
  694.       case OP_GRAHAM_HULL: if (chExists) {
  695.         PSnewpath();
  696.         PSmoveto(X(ch_vertex[0]),Y(ch_vertex[0]));
  697.         for(i=1;i<nch_vertices;i++) PSlineto(X(ch_vertex[i]),Y(ch_vertex[i]));
  698.         PSclosepath();
  699.         PSstroke();
  700.         NXPing();
  701.         break;
  702.         };
  703.       case OP_VORONOI_DIAGRAM: if (vorExists) {
  704.         float *foo;
  705.         NX_MALLOC(foo,float,numPoints*2);
  706.         for(i=0;i<numPoints;i++) {
  707.           int nv,j;
  708.           if ((nv=find_vregion(i,foo)) == -1) {
  709.             NXRunAlertPanel("Voronoi","find_vregion (%d) failed",
  710.                             NULL,NULL,NULL,i);
  711.             break;
  712.             };
  713.           PSnewpath();
  714.           PSmoveto(foo[0],foo[1]);
  715.           for(j=1;j<nv;j++) PSlineto(foo[j*2],foo[j*2+1]);
  716.           PSclosepath();
  717.           PSstroke();
  718.           NXPing();
  719.           };
  720.         NX_FREE(foo);
  721.         break;
  722.         };
  723.       case OP_DELAUNAY_TRIANGULATION: if (vorExists) [self drawDT]; break;
  724.       case OP_GABRIEL_GRAPH: if (ggExists)  [self drawGG]; break;
  725.       case OP_RELATIVE_NEIGHBORHOOD_GRAPH: if (rngExists) [self drawRNG]; break;
  726.       };
  727.     };
  728.   NXPing();
  729.   STATUS("idle");
  730.   return self;
  731. }
  732.  
  733. /* draw an edge with given endpoints */
  734. -drawEdge:(float)x1:(float)y1:(float)x2:(float)y2
  735. { /*assumes focus is locked already and attributes (color, width) already set*/
  736.   PSnewpath();
  737.   PSmoveto(x1,y1);
  738.   PSlineto(x2,y2);
  739.   PSstroke();
  740.   return self;
  741. }
  742.  
  743. - doVoronoiDiagram
  744. {
  745.   NXRect r;
  746.   float xmin,xmax,ymin,ymax;
  747.  
  748.   STATUS("Finding Voronoi diagram");
  749.   [self getBounds:&r];
  750.   xmin=-1.01;
  751.   xmax=1.01;
  752.   ymin=-1.01;
  753.   ymax=1.01;
  754.  
  755.   if (load_vsites(numPoints,data,xmin,ymin,xmax,ymax)) {
  756.     NXRunAlertPanel("Voronoi","load_vsites failed",NULL,NULL,NULL);
  757.     vorExists=NO;
  758.     return self;
  759.     };
  760.   vorExists=YES;
  761.   return self;
  762. }
  763.  
  764. static void checkedge(TwoDView *s,int i1,int i2)
  765. {
  766.   int i;
  767.   int ii1=MIN(i1,i2),ii2=MAX(i1,i2);
  768.   for(i=0;i<s->ndt_edge;i++)
  769.     if ((s->dt_edge[i*2]==ii1)&&(s->dt_edge[i*2+1]==ii2)) return;
  770.   s->dt_edge[s->ndt_edge*2]=ii1;
  771.   s->dt_edge[s->ndt_edge*2+1]=ii2;
  772.   s->ndt_edge++;
  773. }
  774.  
  775. - doDelaunayTriangulation
  776. {
  777.   int ndt,i;
  778.   TRI *tris;
  779.   [self doVoronoiDiagram]; /* cheat! */
  780.  
  781.   STATUS("building Delaunay triangulation data structures");
  782.   /* from list of triangles, obtain minimal list of edged (toss dups) */
  783.   if ((ndt=find_dtriangles(&tris)) == -1) {
  784.     NXRunAlertPanel("Delaunay","find_dtriangles() failed",
  785.                     NULL,NULL,NULL);
  786.     return self;
  787.     };
  788.  
  789.   NX_MALLOC(dt_edge,int,2*3*ndt);
  790.   ndt_edge=0;
  791.  
  792.   for(i=0;i<ndt;i++) {
  793.     int p1=tris[i].v1->pid;
  794.     int p2=tris[i].v2->pid;
  795.     int p3=tris[i].v3->pid;
  796.     checkedge(self,p1,p2);
  797.     checkedge(self,p1,p3);
  798.     checkedge(self,p2,p3);
  799.     };
  800.   dtExists=YES;
  801. #ifdef DEBUG
  802.   for(i=0;i<ndt_edge;i++)
  803.     fprintf(stderr,"%d: %d %d\n",i,dt_edge[i*2],dt_edge[i*2+1]);
  804. #endif
  805.   return self;
  806. }
  807.  
  808.  
  809. - doGabrielGraph
  810. {
  811.   int i;
  812.  
  813.   if (!dtExists)[self doDelaunayTriangulation];
  814.   /* we need to see the Delaunay graph while we're animating */
  815.   NXSetColor([fgColorWell color]);
  816.   PSsetinstance(NO);
  817.   [[[self drawDT] window] flushWindow];
  818.   NXSetColor([highColorWell color]);
  819.   PSsetinstance(YES);
  820.  
  821.   ngg_edge=ndt_edge;
  822.   NX_MALLOC(gg_edge,int,2*ngg_edge);
  823.  
  824.   bcopy(dt_edge,gg_edge,2*ngg_edge*sizeof(int));
  825.  
  826.  
  827.   /*for each DT edge, check circle of influence */
  828.   for(i=0;i<ngg_edge;i++) {
  829.     int i1=gg_edge[i*2], i2=gg_edge[i*2+1];
  830.     int j;
  831.     float x1=X(i1),y1=Y(i1);
  832.     float x2=X(i2),y2=Y(i2);
  833.     float cx=(x1+x2)/2,cy=(y1+y2)/2,r=0.5*sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  834.     char buf[80];
  835. #ifdef DEBUG
  836.     fprintf(stderr,"edge %d (%d %d)\n",i,i1,i2);
  837.     fprintf(stderr,"(%f %f) -> (%f %f)\n",x1,y2,x2,y2);
  838.     fprintf(stderr,"center %f %f rad %f\n",cx,cy,r);
  839. #endif
  840.     sprintf(buf,"considering Delaunay edge %d",i);
  841.     STATUS(buf);
  842.     NI;
  843.     PSnewpath();
  844.     PSarc(cx,cy,r,0.0,360.0);
  845.     PSstroke();
  846.     NXPing();
  847.  
  848.     /* check to see if any other vertex in within the circle of influence.
  849.        if so, remove the edge from the GG */
  850.     for(j=0;j<numPoints;j++) {
  851.       if ((j!=i1)&&(j!=i2)) {
  852.         float x=X(j),y=Y(j),d=(x-cx)*(x-cx)+(y-cy)*(y-cy)-r*r;
  853. #ifdef DEBUG
  854.         fprintf(stderr,"checking %d (%f %f) d %g\n",j,x,y,d);
  855. #endif
  856.         if (d<0.0) { /* clobber edge */
  857.           int k;
  858. #ifdef DEBUG
  859.           fprintf(stderr,"die.\n");
  860. #endif
  861.           for(k=i;k<ngg_edge-1;k++) {
  862.             gg_edge[k*2]=gg_edge[(k+1)*2];
  863.             gg_edge[k*2+1]=gg_edge[(k+1)*2+1];
  864.             };
  865.           ngg_edge--;
  866.           /* highlight the dead edge */
  867.           PSsetinstance(NO);
  868.           PSnewpath();
  869.           PSmoveto(x1,y1);
  870.           PSlineto(x2,y2);
  871.           PSstroke();
  872.           PSsetinstance(YES);
  873.           i--;
  874.           break;
  875.           };
  876.         };
  877.       };
  878.     };
  879.   NX_REALLOC(gg_edge,int,2*ngg_edge);
  880.   ggExists=YES;
  881.   return self;
  882. }
  883.  
  884. - doRelativeNeighborhoodGraph
  885. {
  886.   int i;
  887.  
  888.   if (!dtExists)[self doDelaunayTriangulation];
  889.   /* we need to see the Delaunay graph while we're animating */
  890.   PSsetinstance(NO);
  891.   NXSetColor([fgColorWell color]);
  892.   [[[self drawDT] window] flushWindow];
  893.   NXSetColor([highColorWell color]);
  894.   PSsetinstance(YES);
  895.  
  896.   nrng_edge=ndt_edge;
  897.   NX_MALLOC(rng_edge,int,2*nrng_edge);
  898.  
  899.   bcopy(dt_edge,rng_edge,2*nrng_edge*sizeof(int));
  900.  
  901.  
  902.   /*for each DT edge, check lune of influence */
  903.   for(i=0;i<nrng_edge;i++) {
  904.     int i1=rng_edge[i*2], i2=rng_edge[i*2+1];
  905.     int j;
  906.     float x1=X(i1),y1=Y(i1);
  907.     float x2=X(i2),y2=Y(i2);
  908.     float r2=(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);
  909.     float r=sqrt(r2);
  910.     char buf[80];
  911. #ifdef DEBUG
  912.     fprintf(stderr,"edge %d (%d %d)\n",i,i1,i2);
  913.     fprintf(stderr,"(%f %f) -> (%f %f)\n",x1,y2,x2,y2);
  914.     fprintf(stderr,"rad^2 %f\n",r2);
  915. #endif
  916.     sprintf(buf,"considering Delaunay edge %d",i);
  917.     STATUS(buf);
  918.     NI;
  919.     PSnewpath();
  920.     PSarc(x1,y1,r,0.0,360.0);
  921.     PSmoveto(x2,y2); /* eliminate annoying connector */
  922.     PSarc(x2,y2,r,0.0,360.0);
  923.     PSstroke();
  924.     NXPing();
  925.  
  926.     /* check to see if any other vertex in within the lune of influence.
  927.        if so, remove the edge from the RNG */
  928.     for(j=0;j<numPoints;j++) {
  929.       if ((j!=i1)&&(j!=i2)) {
  930.         float x=X(j),y=Y(j);
  931.         float d1=(x-x1)*(x-x1)+(y-y1)*(y-y1)-r*r;
  932.         float d2=(x-x2)*(x-x2)+(y-y2)*(y-y2)-r*r;
  933. #ifdef DEBUG
  934.         fprintf(stderr,"checking %d (%f %f) d1 %g d2 %g\n",j,x,y,d1,d2);
  935. #endif
  936.         if ((d1<0.0)&&(d2<0.0)) { /* clobber edge */
  937.           int k;
  938. #ifdef DEBUG
  939.           fprintf(stderr,"die.\n");
  940. #endif
  941.           for(k=i;k<nrng_edge-1;k++) {
  942.             rng_edge[k*2]=rng_edge[(k+1)*2];
  943.             rng_edge[k*2+1]=rng_edge[(k+1)*2+1];
  944.             };
  945.           nrng_edge--;
  946.           /* highlight the dead edge */
  947.           PSsetinstance(NO);
  948.           PSnewpath();
  949.           PSmoveto(x1,y1);
  950.           PSlineto(x2,y2);
  951.           PSstroke();
  952.           PSsetinstance(YES);
  953.           i--;
  954.           break;
  955.           };
  956.         };
  957.       };
  958.     };
  959.   NX_REALLOC(rng_edge,int,2*nrng_edge);
  960.   rngExists=YES;
  961.   return self;
  962. }
  963.  
  964. - saveTIFF:sender
  965. {
  966.   NXRect r;
  967.   id bm=[NXBitmapImageRep alloc];
  968.   id sp=[SavePanel new];
  969.   NXStream *stream;
  970.   int fd;
  971.  
  972.   [self lockFocus];
  973.   [self getBounds:&r];
  974.   [bm initData:NULL fromRect:&r];
  975.   [self unlockFocus];
  976.   
  977.   [sp setRequiredFileType:"tiff"];
  978.   if (![sp runModal]) return self;
  979.   fd=open([sp filename],O_CREAT|O_RDWR,0644);
  980.   /*[sp free]; screaming death! screaming death! don't do this! */
  981.   stream = NXOpenFile(fd,NX_READWRITE);
  982.   if (!stream) {
  983.     NXRunAlertPanel("TIFF","NXOpenFile failed",NULL,NULL,NULL);
  984.     perror("NXOpenFile");
  985.     [bm free];
  986.     close(fd);
  987.     return self;
  988.     };
  989.   STATUS("saving image as TIFF");
  990.   [bm writeTIFF:stream usingCompression:NX_TIFF_COMPRESSION_LZW];
  991.   [bm free];
  992.   NXClose(stream);
  993.   close(fd);
  994.   STATUS("idle");
  995.   return self;
  996. }
  997.  
  998. - saveData:sender
  999. {
  1000.   id sp=[SavePanel new];
  1001.   FILE *fp;
  1002.   int i;
  1003.  
  1004.   [sp setRequiredFileType:"2Ddata"];
  1005.   if (![sp runModal]) return self;
  1006.   fp=fopen([sp filename],"w");
  1007.   /*[sp free];*/
  1008.   if (!fp) {
  1009.     NXRunAlertPanel("Data","Couldn't open file",NULL,NULL,NULL);
  1010.     return self;
  1011.     };
  1012.   STATUS("saving data");
  1013.   fprintf(fp,"%d\n",numPoints);
  1014.   for(i=0;i<numPoints;i++) fprintf(fp,"%g %g\n",X(i),Y(i));
  1015.   fclose(fp);
  1016.   STATUS("idle");
  1017.   return self;
  1018. }
  1019.  
  1020. - openData:sender
  1021. {
  1022.   id op=[OpenPanel new];
  1023.   const char *const types[] = { "2Ddata" , NULL};
  1024.   FILE *fp;
  1025.   int i;
  1026.  
  1027.   [op setRequiredFileType:"2Ddata"];
  1028.   [op runModalForDirectory:"." file:"" types:types];
  1029.   chdir([op directory]);
  1030.   if (![op filenames]) {
  1031.     /*[op free];*/
  1032.     return self;
  1033.     };
  1034.  
  1035.   fp=fopen(*[op filenames],"r");
  1036.   /*[op free];*/
  1037.   if (!fp) {
  1038.     NXRunAlertPanel("Open","Couldn't open file",NULL,NULL,NULL);
  1039.     return self;
  1040.     };
  1041.  
  1042.   STATUS("reading data");
  1043.   fscanf(fp,"%d ",&numPoints);
  1044.   initGeometryStuff(self,numPoints);
  1045.   for(i=0;i<numPoints;i++) {
  1046.     fscanf(fp,"%f %f ",&X(i),&Y(i));
  1047.     };
  1048.   fclose(fp);
  1049.   [self disp:0];
  1050.   STATUS("idle");
  1051.   return self;
  1052. }
  1053.  
  1054. - close:sender
  1055. {
  1056.   initGeometryStuff(self,0);
  1057.   [self erase];
  1058.   return self;
  1059. }
  1060.   
  1061. -lineWidthChanged:sender
  1062. {
  1063.   char buf[80];
  1064.   sprintf(buf,"Line width: %6.3f\n",[sender floatValue]);
  1065.   [lineWidthText setStringValue:buf];
  1066.   /* could redraw self here */
  1067.   if (autoUpdate) [self display];
  1068.   return self;
  1069. }
  1070.  
  1071.  
  1072. -setFgColorWell:sender
  1073. {
  1074.   fgColorWell = sender;
  1075.   [sender setColor:NXConvertGrayToColor(0.0)];
  1076.   return self;
  1077. }
  1078.  
  1079. -setBgColorWell:sender
  1080. {
  1081.   bgColorWell = sender;
  1082.   [sender setColor:NXConvertGrayToColor(1.0)];
  1083.   return self;
  1084. }
  1085.  
  1086.  
  1087. -setHighColorWell:sender
  1088. {
  1089.   highColorWell = sender;
  1090.   [sender setColor:NXConvertGrayToColor(0.666)];
  1091.   return self;
  1092. }
  1093.  
  1094. - (BOOL)acceptsFirstResponder { return YES; }
  1095. - (BOOL)acceptsFirstMouse { return YES; }
  1096.  
  1097. - copy:sender
  1098. {
  1099.   id pb=[Pasteboard new];
  1100.   NXStream *stream;
  1101.   char *d;
  1102.   int length,maxLength;
  1103.  
  1104.   STATUS("copying image to clipboard");
  1105.   [pb declareTypes:&NXPostScriptPboard num:1 owner:self];
  1106.   stream=NXOpenMemory(NULL,0,NX_WRITEONLY);
  1107.   [self copyPSCodeInside:NULL to:stream];
  1108.   NXGetMemoryBuffer(stream,&d,&length,&maxLength);
  1109.   [pb writeType:NXPostScriptPboard data:d length:length];
  1110.   NXCloseMemory(stream,NX_FREEBUFFER);
  1111.   [pb free];
  1112.   STATUS("idle");
  1113.   return self;
  1114. }
  1115.  
  1116. - saveEPS:sender
  1117. {
  1118.   NXStream *stream;
  1119.   id sp = [SavePanel new];
  1120.   int fd;
  1121.  
  1122.   [sp setRequiredFileType:"eps"];
  1123.   if (![sp runModal]) return self;
  1124.   fprintf(stderr,"filename %s\n",[sp filename]);
  1125.   fd=open([sp filename],O_CREAT|O_WRONLY,0644);
  1126.   /*[sp free];*/
  1127.   stream = NXOpenFile(fd,NX_WRITEONLY);
  1128.   if (!stream) {
  1129.     NXRunAlertPanel("EPS","NXOpenFile failed",NULL,NULL,NULL);
  1130.     perror("NXOpenFile");
  1131.     return self;
  1132.     };
  1133.   STATUS("saving image as EPS");
  1134.   [self copyPSCodeInside:NULL to:stream];
  1135.   NXClose(stream);
  1136.   close(fd);
  1137.   STATUS("idle");
  1138.   return self;
  1139. }
  1140.  
  1141. @end
  1142.