set here [file dirname [file normalize [info script]]]
::cthulhu::add_directory $here {
}

/*
** This file implements an in-memory balanced binary tree.  This code
** was originally brought in to implement an efficient priority queue
** for the Dijkstra's shortest path algorithm in crewroute.c.
**
** This file contains code imported into READI from another project.
** The text of the original header comment follows:
*/
/*
** A package of routines for handling balanced binary trees.
** Copyright (c) 1990 by D. Richard Hipp
*/

#include <stdlib.h>
#include "odieInt.h"


/* Each node in the binary tree is represented by a single instance
** of the following structure
*/
typedef struct TreeElem TreeElem;
struct TreeElem {
  void *data;             /* Pointer to the user's data */
  void *key;              /* The key associated with this element */
  TreeElem *left;         /* Left daughter */
  TreeElem *right;        /* Right daughter */
  int weight;             /* Weight of this node */
};

/* Turn bulk memory into a Tree structure
*/
void TreeInit(
  Tree *tree,                                 /* Tree object to initialize */
  int (*xCompare)(const void*, const void*),  /* Comparison function */
  void *(*xCopy)(const void*),                /* Key copy function or NULL */
  void (*xFree)(void*)                        /* Key delete function or NULL */
){
  tree->xCompare = xCompare;
  tree->xCopy = xCopy;
  tree->xFree = xFree;
  tree->top = 0;
}

/* Return the number of elements in the tree.
*/
int TreeCount(Tree *pTree){
  if( pTree && pTree->top ){
    return pTree->top->weight;
  }else{
    return 0;
  }
}

/* Delete a single node of the binary tree and all of its children
*/
static void TreeClearNode(TreeElem *p, void (*xFree)(void*)){
  if( p==0 ) return;
  if( p->left ) TreeClearNode(p->left, xFree);
  if( p->right ) TreeClearNode(p->right, xFree);
  if( xFree ){
    xFree(p->key);
  }
  Odie_Free((char *)p);
}

/* Remove all nodes from a tree
*/
void TreeClear(Tree *tree){
  if( tree->top ){
    TreeClearNode(tree->top, tree->xFree);
  }
  tree->top = 0;
}

/* Find the element of the tree (if any) whose key matches "key".
** Return a pointer to the data for that element.  If no match
** is found, return NULL.
*/
void *TreeFind(Tree *tree, const void *key){
  TreeElem *p;

  p = tree->top;
  while( p ){
    int c = tree->xCompare(p->key, key);
    if( c==0 ){
      return p->data;
    }else if( c<0 ){
      p = p->right;
    }else{
      p = p->left;
    }
  }
  return 0;
}

/* If the node with key "key" is the left-most element of the tree, 
** return 0.  If it is the second to the left, return 1.  And so
** forth.
**
** If there is no node in the tree with the key "key", then return
** the number that would have been returned if such a node were
** inserted.
*/
int TreeRank(Tree *tree, const void *key){
  TreeElem *p;
  int rank = 0;

  p = tree->top;
  while( p ){
    int c = tree->xCompare(p->key, key);
    if( c==0 ){
      rank += p->left ? p->left->weight: 0;
      break;
    }else if( c<0 ){
      rank += (p->left ? p->left->weight: 0) + 1;
      p = p->right;
    }else{
      p = p->left;
    }
  }
  return rank;
}

/* Return a pointer to the N-th element of a tree.  (The left-most
** element is number 0, the next is number 1 and so forth.)
*/
static TreeElem *TreeFindNthElem(Tree *tree, int n){
  TreeElem *p;

  p = tree->top;
  while( p ){
    int c = p->left ? p->left->weight : 0;
    if( n==c ){
      return p;
    }else if( n>c ){
      n -= c+1;
      p = p->right;
    }else{
      p = p->left;
    }
  }
  return 0;
}

/* Return the data associated with the N-th element of the tree.  Return
** NULL if there is no N-th element.
*/
void *TreeData(Tree *tree, int n){
  TreeElem *p = TreeFindNthElem(tree,n);
  return p ? p->data : 0;
}

/* Return the key associated with the N-th element of the tree.  Return
** NULL if there is no N-th element.
*/
const void *TreeKey(Tree *tree, int n){
  TreeElem *p = TreeFindNthElem(tree,n);
  if( p ){
    return p->key;
  }else{
    return 0;
  }
}

/*
** Definitions:
** WEIGHT
**   The weight of a node is the total number of children for the node
**   plus 1.  Leaf nodes have a weight of 1.  The root node has a weight
**   which equals the number of nodes in the tree.
**
** BALANCE
**   A node is balanced if the weight of the one child is not more than
**   twice the weight of the other child.
*/

/* The following routine rebalances the tree rooted at *ppElem after
** the insertion or deletion of a single ancestor.
*/
static void TreeBalance(TreeElem **ppElem){
  TreeElem *n;     /* Pointer to self */
  int l,r;         /* Weight of left and right daughters */
  int a,b;         /* Weights of various nodes */

  if( ppElem==0 || (n=*ppElem)==0 ) return;
  l = n->left ? n->left->weight: 0;
  r = n->right ? n->right->weight: 0;
  if( l>r*2 ){    /* Too heavy on the left side */
    TreeElem *nl;     /* Pointer to left daughter */
    TreeElem *nr;     /* Pointer to right daughter */
    int ll, lr;       /* Weight of left daughter's left and right daughter */
    nl = n->left;
    ll = nl->left ? nl->left->weight: 0;
    lr = nl->right ? nl->right->weight: 0;
    if( ll>lr || nl->right==0 ){
      /*
      ** Convert from:  n     to:  nl
      **               / \        /  \
      **              nl  c      a    n
      **             /  \            / \
      **            a    b          b   c
      */
      n->left = nl->right;
      nl->right = n;
      n->weight = a = r + lr + 1;
      nl->weight = a + ll + 1;
      *ppElem = nl;
    }else{
      /*
      ** Convert from:  n        to:  nr
      **               / \           /  \
      **             nl   d        nl    n
      **            /  \          / \   / \
      **           a    nr       a   b c   d
      **               /  \
      **              b    c
      */
      int lrl, lrr;    /* Weight of Great-granddaughter nodes */
      nr = nl->right;
      lrl = nr->left ? nr->left->weight: 0;
      lrr = nr->right ? nr->right->weight: 0;
      nl->right = nr->left;
      nr->left = nl;
      n->left = nr->right;
      nr->right = n;
      n->weight = a = lrr + r + 1;
      nl->weight = b = ll + lrl + 1;
      nr->weight = a + b + 1;
      *ppElem = nr;
    }
  }else if( r>l*2 ){/* Too deep on the right side */
    TreeElem *nl;     /* Pointer to left daughter */
    TreeElem *nr;     /* Pointer to right daughter */
    int rl, rr;       /* Weight of right daughter's left and right daughter */
    nr = n->right;
    rl = nr->left ? nr->left->weight: 0;
    rr = nr->right ? nr->right->weight: 0;
    if( rr>rl || nr->left==0 ){
      /*
      ** Convert from:  n         to:  nr
      **               / \            /  \
      **              a   nr         n    c 
      **                 /  \       / \
      **                b    c     a   b
      */
      n->right = nr->left;
      nr->left = n;
      n->weight = a = l + rl + 1;
      nr->weight = a + rr + 1;
      *ppElem = nr;
    }else{
      /*
      ** Convert from:  n         to:  nl
      **               / \            /  \
      **              a   nr         n    nr
      **                 /  \       / \   / \
      **               nl    d     a   b c   d
      **              /  \
      **             b    c
      */
      int rll,rlr;    /* Weights of great-granddaughter nodes */
      nl = nr->left;
      rll = nl->left ? nl->left->weight: 0;
      rlr = nl->right ? nl->right->weight: 0;
      nr->left = nl->right;
      nl->right = nr;
      n->right = nl->left;
      nl->left = n;
      n->weight = a = l + rll + 1;
      nr->weight = b = rr + rlr + 1;
      nl->weight = a + b + 1;
      *ppElem = nl;
    }
  }else{ /* Node is already balanced.  Just recompute its weight. */
    n->weight = l + r + 1;
  }
}

/* This routine either changes the data on an existing node in the tree,
** or inserts a new node.  "key" identifies the node.  If the data on
** an existing node is changed, then the function returns the old data.
** If a new node is created, NULL is returned.
*/
static void *TreeInsertElement(
  Tree *pTree,                /* The root of the tree */
  void *key,                  /* Insert data at this key */
  void *data                  /* Data to be inserted */
){
  TreeElem *n;
  void *old = 0;
  TreeElem **h[100];  /* Sufficient for a tree with up to 4.0E+17 nodes */
  int level = 0;


  h[0] = &pTree->top;
  level = 1;
  n = pTree->top;
  while( n ){
    int c;
    c = pTree->xCompare(key, n->key);
    if( c<0 ){
      h[level++] = &(n->left);
      n = n->left;
    }else if( c>0 ){
      h[level++] = &(n->right);
      n = n->right;
    }else{
      old = n->data;
      n->data = data;                /* Replace data in an existing node */
      break;
    }
  }
  if( n==0 ){                     /* Insert a leaf node */
    level--;
    n = *h[level] = (TreeElem *)Odie_Alloc( sizeof(TreeElem) );
    if( n==0 ){
      return data;
    }
    n->data = data;
    if( pTree->xCopy ){
      n->key = pTree->xCopy(key);
    }else{
      n->key = key;
    }
    n->left = n->right = 0;
    while( level>=0 ) TreeBalance(h[level--]);
  }
  return old;
}

/* Unlink the N-th node of the tree and return a pointer to that
** node.  (The left-most node is 0, the next node to the right is
** 1 and so forth.)
*/
static TreeElem *TreeDeleteNthElem(TreeElem **ppTop, int N){
  TreeElem *p;        /* For walking the tree */
  int level = 0;      /* Depth of the blancing stack */
  TreeElem **h[100];  /* Balance stack.  100 is sufficient for balancing
                      ** a tree with up to 4.0E+17 nodes */

  h[0] = ppTop;
  level = 1;
  p = *ppTop;
  while( p ){
    int w;
    w = (p->left ? p->left->weight: 0);
    if( N>w ){
      h[level++] = &(p->right);
      p = p->right;
      N -= w+1;
    }else if( N<w ){
      h[level++] = &(p->left);
      p = p->left;
    }else{
      break;
    }
  }
  if( p ){
    level--;
    if( p->left==0 ){
      *h[level] = p->right;
      level--;
    }else if( p->right==0 ){
      *h[level] = p->left;
      level--;
    }else{
      TreeElem *x;
      x = TreeDeleteNthElem(&(p->right),0);
      x->right = p->right;
      x->left = p->left;
      *h[level] = x;
    }
    while( level>=0 ) TreeBalance(h[level--]);
  }
  return p;
}

/* Unlink the node of the tree corresponding to key and return a pointer 
** to that node.
*/
static TreeElem *TreeDeleteElem(Tree *tree, const void *key){
  TreeElem *p;        /* For walking the tree */
  int level = 0;      /* Depth of the blancing stack */
  TreeElem **h[100];  /* Balance stack.  100 is sufficient for balancing
                      ** a tree with up to 4.0E+17 nodes */

  h[0] = &tree->top;
  level = 1;
  p = tree->top;
  while( p ){
    int w;
    w = tree->xCompare(p->key, key);
    if( w<0 ){
      h[level++] = &(p->right);
      p = p->right;
    }else if( w>0 ){
      h[level++] = &(p->left);
      p = p->left;
    }else{
      break;
    }
  }
  if( p ){
    level--;
    if( p->left==0 ){
      *h[level] = p->right;
      level--;
    }else if( p->right==0 ){
      *h[level] = p->left;
      level--;
    }else{
      TreeElem *x;
      x = TreeDeleteNthElem(&(p->right),0);
      x->right = p->right;
      x->left = p->left;
      *h[level] = x;
    }
    while( level>=0 ) TreeBalance(h[level--]);
  }
  return p;
}

/* Insert new data into a node of the tree.  The node is identified
** by "key".
**
** If the new data is NULL, then node is deleted.
**
** If the node aready exists, the new data overwrites the old and
** the old data is returned.  If the node doesn't already exist, then
** a new node is created and the function returns NULL.
*/
void *TreeInsert(Tree *tree, void *key, void *data){
  void *old;
  if( data==0 ){
    TreeElem *elem = TreeDeleteElem(tree, key);
    if( elem ){
      if( tree->xFree ){
        tree->xFree(elem->key);
      }
      old = elem->data;
      Odie_Free((char *)elem);
    }else{
      old = 0;
    }
  }else{
    old = TreeInsertElement(tree,key,data);
  }
  return old;
}

/* Change the data on the n-th node of the tree.  The old data
** is returned.
**
** If data==NULL, then the n-th node of the tree is deleted.  (The
** data associated with that node is still returned.)
**
** If the value N is out-of-bounds, then no new node is created.
** Instead, the "data" parameter is returned.
*/
void *TreeChangeNth(Tree *tree, int n, void *data){
  void *old;
  if( data==0 ){
    TreeElem *elem = TreeDeleteNthElem(&tree->top,n);
    if( elem ){
      if( tree->xFree ){
        tree->xFree(elem->key);
      }
      old = elem->data;
      Odie_Free((char *)elem);
    }else{
      old = 0;
    }
  }else{
    TreeElem *elem = TreeFindNthElem(tree,n);
    if( elem ){
      old = elem->data;
      elem->data = data;
    }else{
      old = data;
    }
  }
  return old;
}

int DLLEXPORT Tree_Init(Tcl_Interp *interp){
  #if IRM_MEM_DEBUG
  Tcl_LinkVar(interp, "module_malloc(tree)", (char*)&nMalloc,
        TCL_LINK_INT | TCL_LINK_READ_ONLY);
  #endif
  return TCL_OK;
}

/*
** Original code by D. Richard Hipp, circa 1990.
** Modified for READI, 2005-11-21.
*/

/* A complete binary tree is defined by an instance of the following
** structure
*/
typedef struct Tree Tree;
struct Tree {
  int (*xCompare)(const void*, const void*); /* Comparison function */
  void *(*xCopy)(const void*);               /* Key copy function, or NULL */
  void (*xFree)(void*);                      /* Key delete function */
  struct TreeElem *top;                      /* The top-most node of the tree */
};

void TreeInit(
  Tree *tree,                                 /* Tree object to initialize */
  int (*xCompare)(const void*, const void*),  /* Comparison function */
  void *(*xCopy)(const void*),                /* Key copy function or NULL */
  void (*xFree)(void*)                        /* Key delete function or NULL */
);
void TreeClear(Tree *tree);
void *TreeChangeNth(Tree *tree, int n, void *data);
void *TreeInsert(Tree *tree, void *key, void *data);
void *TreeFind(Tree *tree, const void *key);
int TreeRank(Tree *tree, const void *key);
void *TreeData(Tree *tree, int n);
const void *TreeKey(Tree *tree, int n);
int TreeCount(Tree*);

set here [file dirname [file normalize [info script]]]
foreach file {plotter.c slicer.c wallset.c} {
  ::cthulhu::detect_cases [file join $here generic $file]
}
::cthulhu::add_directory [file join $here generic] {}

/*
** Load the entire geometry module
*/
#include "odieInt.h"

/*
** Declare our Tcl Obj Types
*/

const Tcl_ObjType polygon_tclobjtype = {
  "polygon", /* name */
  &PolygonObj_freeIntRepProc, /* freeIntRepProc */
  &PolygonObj_dupIntRepProc, /* dupIntRepProc */
  &PolygonObj_updateStringProc, /* updateStringProc */
  &PolygonObj_setFromAnyProc /* setFromAnyProc */
};

const Tcl_ObjType segmentset_tclobjtype = {
  "segmentset", /* name */
  &SegmentSetObj_freeIntRepProc, /* freeIntRepProc */
  &SegmentSetObj_dupIntRepProc, /* dupIntRepProc */
  &SegmentSetObj_updateStringProc, /* updateStringProc */
  &SegmentSetObj_setFromAnyProc /* setFromAnyProc */
};

Tcl_Obj *Odie_NewPolygonObj(Poly *pPoly) {
  Tcl_Obj *pResult;
  pResult=Tcl_NewObj();
  Tcl_InvalidateStringRep(pResult);
  pResult->internalRep.otherValuePtr=pPoly;
  pResult->typePtr=&polygon_tclobjtype;
  return pResult;
}

int Odie_GetPolygonFromObj(Tcl_Interp *interp,Tcl_Obj *objPtr,Poly **ptr,int *created) {
  Poly *p;
  *created=0;

  if(objPtr->typePtr) {
    if(objPtr->typePtr==&polygon_tclobjtype && objPtr->internalRep.otherValuePtr) {
      /*
      ** Object is already of the type requested
      */
      *ptr=objPtr->internalRep.otherValuePtr;
      return TCL_OK;
    }
  }
  int k,i;
  if( Tcl_ListObjLength(interp, objPtr, &k) ) return TCL_ERROR;
  if( k<6 ){
    Tcl_AppendResult(interp, "need at least 3 vertices", 0);
    return TCL_ERROR;
  }
  if( k&1 ){
    Tcl_AppendResult(interp, "coordinates should come in pairs", 0);
    return TCL_ERROR;
  }
  p=(Poly *)Odie_Alloc(sizeof(*p)+(k+2)*sizeof(p->v[0]));
  p->nVertex=k/2;
  for(i=0; i<p->nVertex; i++){
    Tcl_Obj *pElem;
    double d;
    Tcl_ListObjIndex(0, objPtr, i*2, &pElem);
    if(Tcl_GetDoubleFromObj(interp, pElem, &d)) goto createfail;
    p->v[i][X_IDX] = d;
    Tcl_ListObjIndex(0, objPtr, i*2+1, &pElem);
    if(Tcl_GetDoubleFromObj(interp, pElem, &d)) goto createfail;
    p->v[i][Y_IDX] = d;
  }

  if(Odie_PolygonComputeArea(interp,p)==TCL_OK) {
    *ptr=p;
    if(Tcl_IsShared(objPtr)) {
      *created=1;
    } else {
      /* Shimmer this object to the requested type */
      if(objPtr->typePtr && objPtr->typePtr->freeIntRepProc) {
        Tcl_FreeInternalRepProc *freeIntRepProc=objPtr->typePtr->freeIntRepProc;
        freeIntRepProc(objPtr);
      }
      Tcl_InvalidateStringRep(objPtr);
      objPtr->internalRep.otherValuePtr=p;
      objPtr->typePtr=&polygon_tclobjtype;
    }
    return TCL_OK;
  }
  
createfail:
  Odie_Free((char *)p);
  return TCL_ERROR;
}

int PolygonObj_setFromAnyProc(Tcl_Interp *interp,Tcl_Obj *objPtr) {
  if(objPtr->typePtr) {
    if(objPtr->typePtr==&polygon_tclobjtype) {
      /*
      ** Object is already of the type requested
      */
      return TCL_OK;
    }
  }
  Poly *p;
  int created=0;
  if(Odie_GetPolygonFromObj(interp,objPtr,&p,&created)) {
    return TCL_ERROR;
  }
  Tcl_InvalidateStringRep(objPtr);
  objPtr->internalRep.otherValuePtr=p;
  objPtr->typePtr=&polygon_tclobjtype;
  return TCL_OK;
}

void PolygonObj_updateStringProc(Tcl_Obj *objPtr) {
  char outbuffer[128];
  Tcl_DString result;
  Poly *p=objPtr->internalRep.otherValuePtr;
  int i,j;
  Tcl_DStringInit(&result);
  j=p->nVertex-1;
  if(p->v[0][X_IDX]==p->v[j][X_IDX] && p->v[0][Y_IDX]==p->v[j][Y_IDX]) {
    j=p->nVertex-2;
  }
  for(i=0; i<=j; i++){
    sprintf(outbuffer,"%g %g",(float)p->v[i][X_IDX],(float)p->v[i][Y_IDX]);
    Tcl_DStringAppendElement(&result,outbuffer);
  }
  objPtr->length=Tcl_DStringLength(&result);
  objPtr->bytes=Odie_Alloc(objPtr->length+1);
  strcpy(objPtr->bytes,Tcl_DStringValue(&result));
  Tcl_DStringFree(&result);
}

void PolygonObj_dupIntRepProc(Tcl_Obj *srcPtr,Tcl_Obj *dupPtr) {
  Poly *src=srcPtr->internalRep.otherValuePtr;
  int size=sizeof(*src)+src->nVertex*sizeof(src->v[0]);
  Poly *copy=(Poly *)Odie_Alloc(size);

  memcpy(copy,src,size);
  Tcl_InvalidateStringRep(dupPtr);
  dupPtr->typePtr=&polygon_tclobjtype;
  dupPtr->internalRep.otherValuePtr=copy;
}

void PolygonObj_freeIntRepProc(Tcl_Obj *objPtr) {
  Odie_Free(objPtr->internalRep.otherValuePtr);
  objPtr->internalRep.otherValuePtr=NULL;
  objPtr->typePtr=NULL;
}

void Segset_Insert_Polygon(SegSet *pSet,Poly *p,int fill) {
  int i;
  if(p->nVertex>0) {
    VECTORXY *P;
    P=&p->v[0];
    for(i=1; i<p->nVertex; i++){
      SegSetInsert(pSet,*P,p->v[i],1);
      P=&p->v[i];
    }
    SegSetInsert(pSet,*P,p->v[0],1);
  }
}

void SegSetCopy(SegSet *dest,SegSet *src) {
  SegSetClear(dest);
  memset(dest, 0, sizeof(SegSet));

  Link *pLoop, *pNext;
  for(pLoop=src->pAll; pLoop; pLoop=pNext){
    Segment *pAB;
    pAB = pLoop->pLinkNode;
    pNext = pLoop->pNext;
    SegSetInsert(dest,pAB->from,pAB->to,pAB->isBoundary);
  }
}


/*
** Find and return the line segment that goes from A to B.  Return NULL
** if there is not such line segment
*/
Segment *SegSetFind(SegSet *pSet, VectorXY A, VectorXY B){
  Link *pX;
  Segment *p;
  int h;
  h = hashVectorXY(A);
  for(pX=pSet->hashFrom[h]; pX; pX=pX->pNext){
    p = pX->pLinkNode;
    if( sameVectorXY(p->from, A) && sameVectorXY(p->to, B) ){
      return p;
    }
  }
  return 0;
}

Tcl_Obj *Odie_NewSegmentSetObj(SegSet *pSegSet) {
  Tcl_Obj *pResult;
  pResult=Tcl_NewObj();
  Tcl_InvalidateStringRep(pResult);
  pResult->internalRep.otherValuePtr=pSegSet;
  pResult->typePtr=&segmentset_tclobjtype;
  return pResult;
}

int Odie_GetSegmentSetFromObj(
  Tcl_Interp *interp,
  Tcl_Obj *objPtr,
  SegSet **ptr,
  int *created
) {
  SegSet *pSet=NULL;

  *created=0;  
  if(objPtr->typePtr && objPtr->typePtr->setFromAnyProc==&SegmentSetObj_setFromAnyProc) {
    /*
    ** Object is already of the type requested
    */
    *ptr=objPtr->internalRep.otherValuePtr;
    return TCL_OK;
  }
  if(objPtr->typePtr && objPtr->typePtr->setFromAnyProc==&PolygonObj_setFromAnyProc) {
    *created=1;
    /*
    ** Convert from a polygon
    */
    pSet=(SegSet *)Odie_Alloc(sizeof(SegSet));
    Poly *p=objPtr->internalRep.otherValuePtr;
    Segset_Insert_Polygon(pSet,p,1);
    *ptr=pSet;
    return TCL_OK;
  }

  int i,n;
  if( Tcl_ListObjLength(interp, objPtr, &n) ) return TCL_ERROR;
  if( n%4!=0 ){
    Tcl_AppendResult(interp, "VECTORS argument should contain a multiple of 4 values", 0);
    return TCL_ERROR;
  }
  pSet=(SegSet *)Odie_Alloc(sizeof(SegSet));

  for(i=0; i<n; i+=4){
    int j;
    Segment *found;
    double x[4];
    VECTORXY A,B;

    for(j=0; j<4; j++){
      Tcl_Obj *pObj;
      Tcl_ListObjIndex(0, objPtr, i+j, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &x[j]) ){
        goto createfail;
      }
    }
    A[X_IDX]=x[0];
    A[Y_IDX]=x[1];
    B[X_IDX]=x[2];
    B[Y_IDX]=x[3];
    /*
    ** Do not insert a vector into the wallset if it
    ** matches a vector already given. It's either redundent
    ** or the edge of a hole
    */
    found=SegSetFind(pSet,A,B);
    if(found) continue;
    SegSetInsert(pSet, A, B, 1);
  }
  
  //Link *pLoop, *pNext;
  //for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
  //  Segment *pAB;
  //  pAB = pLoop->pLinkNode;
  //  pNext = pLoop->pNext;
  //  int h = hashPoint(pAB->from);
  //}
  
  *created=1;
  *ptr=pSet;
  return TCL_OK;
  
createfail:
  SegSetClear(pSet);

  Odie_Free((char *)pSet);
  return TCL_ERROR;
}

int SegmentSetObj_setFromAnyProc(Tcl_Interp *interp,Tcl_Obj *objPtr) {
  if(objPtr->typePtr) {
    if(objPtr->typePtr==&segmentset_tclobjtype) {
      /*
      ** Object is already of the type requested
      */
      return TCL_OK;
    }
  }
  SegSet *p;
  int created=0;
  if(Odie_GetSegmentSetFromObj(interp,objPtr,&p,&created)) {
    return TCL_ERROR;
  }
  Tcl_InvalidateStringRep(objPtr);
  objPtr->internalRep.otherValuePtr=p;
  objPtr->typePtr=&segmentset_tclobjtype;
  return TCL_OK;
}

void SegmentSetObj_updateStringProc(Tcl_Obj *objPtr) {
  char outbuffer[128];
  Tcl_DString result;
  SegSet *pSet=objPtr->internalRep.otherValuePtr;
  Tcl_DStringInit(&result);
  Link *pLoop, *pNext;
  for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
    Segment *pAB;
    pAB = pLoop->pLinkNode;
    pNext = pLoop->pNext;
    sprintf(outbuffer,"%g %g %g %g %d",(float)pAB->from[X_IDX],(float)pAB->from[Y_IDX],(float)pAB->to[X_IDX],(float)pAB->to[Y_IDX],pAB->isBoundary);
    Tcl_DStringAppendElement(&result,outbuffer);
  }
  objPtr->length=Tcl_DStringLength(&result);
  objPtr->bytes=Odie_Alloc(objPtr->length+1);
  strcpy(objPtr->bytes,Tcl_DStringValue(&result));
  Tcl_DStringFree(&result);
}

void SegmentSetObj_dupIntRepProc(Tcl_Obj *srcPtr,Tcl_Obj *dupPtr) {
  SegSet *src=(SegSet*)srcPtr->internalRep.otherValuePtr;
  SegSet *dest=(SegSet*)Odie_Alloc(sizeof(SegSet));

  SegSetCopy(dest,src);

  dupPtr->typePtr=srcPtr->typePtr;
  dupPtr->internalRep.otherValuePtr=dest;
}

void SegmentSetObj_freeIntRepProc(Tcl_Obj *objPtr) {
  SegSet *set=(SegSet *)objPtr->internalRep.otherValuePtr;
  SegSetClear(set);
  Odie_Free(objPtr->internalRep.otherValuePtr);
}


int Odie_GetSegmentGetFromVar(
  Tcl_Interp *interp,
  Tcl_Obj *varName,
  SegSet **dest
) {
  Tcl_Obj *objPtr=Tcl_ObjGetVar2(interp,varName,NULL,0);
  if(!objPtr) {
    Tcl_ResetResult(interp);
    *dest=(SegSet *)Odie_Alloc(sizeof(SegSet));
    return TCL_OK;
  }
  int created;
  SegSet *src;
  if( Odie_GetSegmentSetFromObj(interp, objPtr, &src, &created) ) return TCL_ERROR;
  *dest=(SegSet *)Odie_Alloc(sizeof(SegSet));
  SegSetCopy(*dest,src);
  return TCL_OK;
}

DLLEXPORT int Odie_Geometry_Init(Tcl_Interp *interp) {

  Tcl_RegisterObjType(&polygon_tclobjtype);
  Tcl_RegisterObjType(&segmentset_tclobjtype);
  
  if(Odie_Segset_Init(interp)) return TCL_ERROR;
  if(Odie_Polygon_Init(interp)) return TCL_ERROR;
  if(Odie_Shapes_Init(interp)) return TCL_ERROR;
  if(Odie_Mathtools_Init(interp)) return TCL_ERROR;
  
  Tcl_CreateObjCommand(interp, "plotter", Odie_PlotterCreateProc, 0, 0);
  Tcl_CreateObjCommand(interp, "slicer", Odie_SlicerCreateProc, 0, 0);
  Tcl_CreateObjCommand(interp, "wallset", Odie_WallsetCreateProc, 0, 0);
  return TCL_OK;
}

/*
** A bounding box.
*/
typedef unsigned char u8;

/*
** Size of the hash tables
*/
#define SZ_HASH  419

/*
** Grain size for rounding
*/
#define GRAIN 2.5


/*
** An instance of the following structure is used as an entry in
** a doubly linked list.
*/
typedef struct Link Link;
struct Link {
  Link *pNext, **ppPrev; /* Next and prior nodes on the list */
  void *pLinkNode;    /* Structure that this link is a part of */
};

typedef struct Box Box;
struct Box {
  double t, b, l, r;
};

/*
** Every polygon is an instance of the following structure:
*/
typedef struct Point Vertex;
typedef struct Point Point;

struct Point { double x, y; };

typedef struct Poly Poly;
struct Poly {
  int id;            /* Numeric ID of the polygon */
  int nVertex;        /* Number of verticies */
  double area;        /* Area contained within the polygon */
  Box bbox;           /* The bounding box */
  double v[][2];        /* Set of vertices */          
};

/*
** Each straight line in the wallset is recorded as an instance in the
** following structure.
**
** Coordinates are always stored as integers in a right-handed coordinate
** system.  Multiply external coordinates by Wallset.rXZoom and
** Wallset.rYZoom to get the internal coordinates stored here.  Divide
** the coordinates stored here by the zoom factors to get external
** coordinates.
*/

typedef struct Segment Segment;
struct Segment {
  int id;
  double from[2];       /* Beginning coordinate */
  double to[2];        /* Ending coordinate */
  int idLC, idRC;   /* ID numbers of compartments to the left and right */

  Link pAll;         /* All segments */
  Link pHash;        /* All segments with the same hash on id */
  Link pFrom;        /* All segments with the same hash on from */
  Link pTo;          /* All segments with the same hash on to */
  Link pSet;         /* A temporary subset of segments */

  double score;     /* Candidate vertex score */
  unsigned int notOblique:1; /* True if next segment is not at an oblique angle */
  unsigned int isBoundary:8; /* True if this is a boundary segment */
  unsigned int ignore:1;        /* Causes this segment to be ignored on some operations */
  unsigned int midpoint:1;        /* Causes this segment to be ignored on some operations */
  int isRight:4; /* -1 isleft - 0 straight - 1 right */
};

/*
** A complete set of segments
*/
typedef struct SegSet SegSet;
struct SegSet {
  int shared;
  int nSeg;              /* Number of segments in the set */
  Segment *pCurrent;     /* Current segment */
  Link *pAll;            /* All segments connected by Segment.all */
  Link *hashFrom[SZ_HASH];  /* Hash on Segment.orig */
};

/*
** A boundary consists of three or more segments.  Each segment contains
** a direction flag.  The following structure holds a single element
** of a boundary.
*/
typedef struct Boundary Boundary;
struct Boundary {
  Segment *pSeg;   /* The segment of the boundary */
  int backwards;   /* True if the boundary moves backwards on the segment */
};

/*
** Instances of the following structure are used to speed up the
** "WS primary" method.  Each instance of this structure records
** a primary wall for a compartment and a rectangular bounding box
** for that compartment.  When searching for a compartment that
** contains a point, we can do a quick scan of the bounding boxes
** in linear time.
*/
typedef struct ComptBox ComptBox;
struct ComptBox {
  ComptBox *pNext;     /* Next ComptBox in a list of them all */
  Box bbox;            /* The bounding box */
  double area;         /* Area of the bounding box, used for sorting */
  Boundary prim;       /* Primary boundary wall for the compartment */
};

/*
** A wallset is an instance of this structure
*/
typedef struct Wallset Wallset;
struct Wallset {
  int busy;                 /* Do not delete or insert if true */
  double rXZoom, rYZoom;    /* Zoom for input and output */
  ComptBox *pComptBox;      /* Cache of compartment boxes */
  Link *pAll;               /* Any of the segments in the Segment.all ring */
  Link *hashId[SZ_HASH];    /* Hash table for Segment.id */
  Link *hashFrom[SZ_HASH];  /* Hash table for Segment.x0,Segment.y0 */
  Link *hashTo[SZ_HASH];    /* Hash table for Segment.x1,Segment.y1 */
};

/*
** A slicer is an instance of the following structure.
*/
typedef struct Slicer Slicer;
struct Slicer {
  int nSlice;        /* Number of slices */
  struct OneSlice {  /* One entry per slice */
    int idx;            /* Index of this slice in a[] */
    int did;            /* Integer deck id */
    int above;            /* Integer deck id above */
    int below;            /* Integer deck id below */
    
    char *zName;        /* Name of this slice */
    double z;           /* Z coordinate of this slice */
    int nXZ;            /* Number of entries in xz[] */
    double rXShift;     /* Change X values by this amount times rZoom */
    double *xz;         /* Alternating X and Z profile values */
    double mnX, mxX;    /* Min and max Y for the slice (actual coord space) */
    double mnY, mxY;    /* Min and max Y for the slice (actual coord space) */
    double top, btm;    /* Top and bottom of slice in canvas coords. top<btm */
    double upperbound;  /* Half-way from a[i].top to a[i+1].btm */
  } *a;              /* Entries malloced and sorted by increasing z */
  double upper_height;  /* Height to assume for top level if given negative coords */
  double rZoom;      /* Multiply canvas coord by this to get actual coord */
  double rXOffset;   /* X-Shift amount per unit of height */
};

/*
** Initialize a new list element so that it is on a list all by itself
** and so that the ring is a member of structure pContainer.
*/
#define LinkInit(R,C)  ((R).pNext = 0, (R).ppPrev = 0, (R).pLinkNode = (void *)(C))

/*
** Routines for handling a doubly linked list
*/
#include "odieInt.h"

/*
** Remove a link element from its list.
*/
CTHULHU_INLINE void LinkRemove(Link *p){
  assert( p->pLinkNode !=0 );
  if( p->ppPrev ){
    assert( *p->ppPrev==p );
    *p->ppPrev = p->pNext;
  }
  if( p->pNext ){
    assert( p->pNext->ppPrev == &p->pNext );
    p->pNext->ppPrev = p->ppPrev;
  }
  p->pNext = 0;
  p->ppPrev = 0;
}

/* 
** Add a link to a list
*/
CTHULHU_INLINE void LinkInsert(Link **ppRoot, Link *pLink){
  assert( pLink->ppPrev==0 );
  assert( pLink->pNext==0 );
  assert( pLink->pLinkNode!=0 );
  pLink->ppPrev = ppRoot;
  pLink->pNext = *ppRoot;
  if( pLink->pNext ){
    pLink->pNext->ppPrev = &pLink->pNext;
  }
  *ppRoot = pLink;
}

/*
** Return the number of elements on a linked list
*/

CTHULHU_INLINE int LinkCount(Link *pRoot){
  int cnt = 0;
  while( pRoot ){
    cnt++;
    pRoot = pRoot->pNext;
  }
  return cnt;
}

/*
** Compute a hash on an integer.
*/
CTHULHU_INLINE int hashInt(int x){
  return (x & 0x7fffffff) % SZ_HASH;
}

/*
** Round a coordinate to its nearest grain boundary
*/
CTHULHU_INLINE long intCoord(double x) {
  long idxX = x/GRAIN + (x>0.0 ? 0.5 : -0.5);
  return idxX*GRAIN;
}

/*
** Compute a hash on a point.
*/
CTHULHU_INLINE int hashPoint(VECTORXY p){
  int idxX = p[X_IDX]/GRAIN;
  int idxY = p[Y_IDX]/GRAIN;
  return hashInt(idxX+idxY);
}

CTHULHU_INLINE int hashVectorXY(VECTORXY p){
  int idxX = p[X_IDX]/GRAIN;
  int idxY = p[Y_IDX]/GRAIN;
  return hashInt(idxX+idxY);
}


CTHULHU_INLINE double roundCoord(double x){
    return intCoord(x);
}

/*
** Compute a hash on a pair of floating point number.
*/
CTHULHU_INLINE int hashCoord(double x, double y){
  long idxX = intCoord(x);
  long idxY = intCoord(y);
  return hashInt(idxX+idxY);
}

/*
** Compare to floating point values and return negative, zero, or
** positive if the first value is less than, equal to, or greater
** than the second.
*/
CTHULHU_INLINE int floatCompare(double x0, double x1){
  double diff = x1 - x0;
  if( diff>-GRAIN && diff<GRAIN ){
    diff = 0.0;
  }
  return (int)diff;
}

/*
** Return TRUE if x0,y0 is the same point as x1,y1
*/
CTHULHU_INLINE int samePoint(double x0, double y0, double x1, double y1){
  return floatCompare(x0,x1)==0 && floatCompare(y0,y1)==0;
}

CTHULHU_INLINE int sameVectorXY(VectorXY A, VectorXY B){
  return floatCompare(A[X_IDX],B[X_IDX])==0 && floatCompare(A[Y_IDX],B[Y_IDX])==0;
}

/*
** This file is machine generated. Changes will
** be overwritten on the next run of cstruct.tcl
*/
#include "odieInt.h"
#define ODIE_REAL_TOLERANCE 0.001

CTHULHU_INLINE int ODIE_Real_Is_Zero(double x) {
  if(fabs(x) < __FLT_EPSILON__) {
    return 1;
  }
  return 0;
}

int ODIE_Fuzzy_Compare_TclObj(Tcl_Obj *avalue,Tcl_Obj *bvalue) {
  int result;
  double ax,bx;
  if(Tcl_GetDoubleFromObj(NULL,avalue,&ax)) {
    char *a,*b;
    a=Tcl_GetString(avalue);
    b=Tcl_GetString(bvalue);
    result=strcmp(a,b);
  } else if(Tcl_GetDoubleFromObj(NULL,bvalue,&bx)) {
    char *a,*b;
    a=Tcl_GetString(avalue);
    b=Tcl_GetString(bvalue);
    result=strcmp(a,b);
  } else {
    result=ODIE_Fuzzy_Compare(ax,bx);
  }
  return result;
}

Tcl_Obj *ODIE_NewFuzzyObj(double ax) {
  double bx;
  const char *use_format="%g";
  char newstring[128];
  
  if(ODIE_Real_Is_Zero(ax)) {
    return ODIE_REAL_ZERO();
  }
  bx=abs(ax);
  if(bx<1e-4 || bx>1e5) {
    use_format="%g";
  } else if(bx>10.0) {
    use_format="%.1f";
  } else if (bx>1.0) {
    use_format="%.2f";
  } else {
    use_format="%.3f";
  }
  sprintf(newstring,use_format,ax);
  return Tcl_NewStringObj(newstring,-1);
}

CTHULHU_INLINE int ODIE_Fuzzy_Compare(double avalue,double bvalue) {
  /* Handle the simple cases */
  double c=bvalue-avalue;
  if(fabs(c) < __FLT_EPSILON__) return 0;
  if(avalue>bvalue) return 1;
  return -1;
}

CTHULHU_INLINE int ODIE_Fuzzy_GTE(double avalue,double bvalue) {
  /* Handle the simple cases */
  if (avalue==bvalue) return 1;
  if (avalue<=ODIE_REAL_TOLERANCE && bvalue>ODIE_REAL_TOLERANCE) return 0;
  if (avalue>bvalue) return 1;

  /* Add epsilon to the send */
  avalue+=ODIE_REAL_TOLERANCE;
  if (avalue>=bvalue) return 2;

  /* For large quantities, loose the decimal points 
  if(avalue>100.0 && bvalue>100.0) {
    avalue=ceil(avalue);
    bvalue=floor(bvalue);
    if (avalue>=bvalue) return 2;
  }
  */
  return 0;
}

/* Detect of two lines are colinear */
CTHULHU_INLINE int Odie_IsColinear(double x1,double y1,double x2,double y2,double x3,double y3) {
  double c=(x3-x1)*(y2-y1)-(y3-y1)*(x2-x1);
  if(fabs(c) < __FLT_EPSILON__) return 1;
  return 0;
}

/*
** Detect the intersection of two lines
** Returns:
** 0 - no overlap
** 1 - AX1 is on line BX1-BY1
** 2 - AX2 is on line BX1-BY1
** 4 - BX1 is on line AX1-AX2
** 8 - BX2 is on line AX1-AX2
*/
CTHULHU_INLINE int ODIE_Math_LineLineCoincident(
double ax1, double ay1,
double ax2, double ay2,
double bx1, double by1,
double bx2, double by2)
{
  double denom,numera,numerb;

  denom  = (by2-by1) * (ax2-ax1) - (bx2-bx1) * (ay2-ay1);   
  /* Are the line parallel */
  if (!ODIE_Real_Is_Zero(denom)) {
   return 0;
  }
  numera = (bx2-bx1) * (ay1-by1) - (by2-by1) * (ax1-bx1);
  numerb = (ax2-ax1) * (ay1-by1) - (ay2-ay1) * (ax1-bx1);
   
  if (!ODIE_Real_Is_Zero(numera) || !ODIE_Real_Is_Zero(numerb)) {
    return 0;
  }
  return 1;
#ifdef NEVER
  VectorXY A,B,C,D;
  A[X_IDX]=ax1;
  A[Y_IDX]=ay1;
  B[X_IDX]=ax2;
  B[Y_IDX]=ay2;
  C[X_IDX]=bx1;
  C[Y_IDX]=by1;
  D[X_IDX]=bx2;
  D[Y_IDX]=by2;

  int result=0;
  if(ODIE_Math_PointOnSegment(C,D,A)) {
    result|=1;
  }
  if(ODIE_Math_PointOnSegment(C,D,B)) {
    result|=2;
  }
  if(ODIE_Math_PointOnSegment(A,B,C)) {
    result|=4;
  }
  if(ODIE_Math_PointOnSegment(A,B,D)) {
    result|=8;
  }
  return result;
#endif
}

/*
** Detect the intersection of two lines
** Adapted from: http://paulbourke.net/geometry/lineline2d/pdb.c
*/
int ODIE_Math_LineLineIntersect(
double ax1, double ay1,
double ax2, double ay2,
double bx1, double by1,
double bx2, double by2,
double *x, double *y)
{
   double mua,mub;
   double denom,numera,numerb;


   denom  = (by2-by1) * (ax2-ax1) - (bx2-bx1) * (ay2-ay1);
   numera = (bx2-bx1) * (ay1-by1) - (by2-by1) * (ax1-bx1);
   numerb = (ax2-ax1) * (ay1-by1) - (ay2-ay1) * (ax1-bx1);
   
   /* Are the line parallel */
   if (ODIE_Real_Is_Zero(denom)) {
    if (ODIE_Real_Is_Zero(numera) && ODIE_Real_Is_Zero(numerb)) {
      /* Are the line coincident? */
      int within=1;
      if(ax2>ax1) {
        if(bx1>ax2 && bx2>ax2) {
          within=0;
        } else if(bx1<ax1 && bx2<ax1) {
          within=0;
        }
      } else {
        if(bx1>ax1 && bx2>ax1) {
          within=0;
        } else if(bx1<ax2 && bx2<ax2) {
          within=0;
        }
      }
      if(ay2>ay1) {
        if(by1>ay2 && by2>ay2) {
          within=0;
        } else if(by1<ay1 && by2<ay1) {
          within=0;
        }
      } else {
        if(by1>ay1 && by2>ay1) {
          within=0;
        } else if(by1<ay2 && by2<ay2) {
          within=0;
        }
      }
      if(within) {
        *x = (ax1 + ax2) / 2;
        *y = (ay1 + ay2) / 2;
        return(1);
      }
    }
    *x = 0;
    *y = 0;
    return(0);
   }

   /* Is the intersection along the the segments */
   mua = numera / denom;
   mub = numerb / denom;
   if (mua < 0 || mua > 1 || mub < 0 || mub > 1) {
    *x = 0;
    *y = 0;
    return(0);
   }
   *x = ax1 + mua * (ax2 - ax1);
   *y = ay1 + mua * (ay2 - ay1);
   return(1);
}

/*
** Detect the intersection of a line and a sphere
** Adapted from: http://http://paulbourke.net/geometry/circlesphere/raysphere.c
*/
int ODIE_Math_LineSphereIntersect(
double p1_x, double p1_y, double p1_z,
double p2_x, double p2_y, double p2_z,
double sc_x, double sc_y, double sc_z,
double r,
double *mu1, double *mu2)
{
   double a,b,c;
   double bb4ac;
   double dp_x,dp_y,dp_z;
   *mu1 = 0;
   *mu2 = 0;

   dp_x = p2_x - p1_x;
   dp_y = p2_y - p1_y;
   dp_z = p2_z - p1_z;
   a = dp_x * dp_x + dp_y * dp_y + dp_z * dp_z;
   b = 2 * (dp_x * (p1_x - sc_x) + dp_y * (p1_y - sc_y) + dp_z * (p1_z - sc_z));
   c = sc_x * sc_x + sc_y * sc_y + sc_z * sc_z;
   c += p1_x * p1_x + p1_y * p1_y + p1_z * p1_z;
   c -= 2 * (sc_x * p1_x + sc_y * p1_y + sc_z * p1_z);
   c -= r * r;
   bb4ac = b * b - 4 * a * c;
   if (ODIE_Real_Is_Zero(a) || bb4ac < 0) {
      return(0);
   }

   *mu1 = (-b + sqrt(bb4ac)) / (2 * a);
   *mu2 = (-b - sqrt(bb4ac)) / (2 * a);

   return(1);
}
/*
** Detect the intersection of a line and a circle
** Adapted from: http://http://paulbourke.net/geometry/circlesphere/raysphere.c
*/
int ODIE_Math_LineCircleIntersect(
double p1_x, double p1_y,
double p2_x, double p2_y,
double sc_x, double sc_y,
double r,
double *mu1, double *mu2)
{
   double a,b,c;
   double bb4ac;
   double dp_x,dp_y;
   *mu1 = 0;
   *mu2 = 0;

   dp_x = p2_x - p1_x;
   dp_y = p2_y - p1_y;
   a = dp_x * dp_x + dp_y * dp_y;
   b = 2 * (dp_x * (p1_x - sc_x) + dp_y * (p1_y - sc_y));
   c = sc_x * sc_x + sc_y * sc_y;
   c += p1_x * p1_x + p1_y * p1_y;
   c -= 2 * (sc_x * p1_x + sc_y * p1_y);
   c -= r * r;
   bb4ac = b * b - 4 * a * c;
   if (ODIE_Real_Is_Zero(a) || bb4ac < 0) {
      return(0);
   }

   *mu1 = (-b + sqrt(bb4ac)) / (2 * a);
   *mu2 = (-b - sqrt(bb4ac)) / (2 * a);

   return(1);
}

/*
** This file is copyright Test and Evaluation Solutions, LLC
** See license.terms for details of usage
*/


static int  odiemath_method_colinear (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  if( objc != 7 ){
      Tcl_WrongNumArgs(interp, 1, objv, "colinear x1 y1 x2 y2 x3 y3");
      return TCL_ERROR;
  }
  double x1,y1,x2,y2,x3,y3;
  if (Tcl_GetDoubleFromObj(interp,objv[1],&x1)) return TCL_ERROR;
  if (Tcl_GetDoubleFromObj(interp,objv[2],&y1)) return TCL_ERROR;
  if (Tcl_GetDoubleFromObj(interp,objv[3],&x2)) return TCL_ERROR;
  if (Tcl_GetDoubleFromObj(interp,objv[4],&y2)) return TCL_ERROR;
  if (Tcl_GetDoubleFromObj(interp,objv[5],&x3)) return TCL_ERROR;
  if (Tcl_GetDoubleFromObj(interp,objv[6],&y3)) return TCL_ERROR;
  Tcl_SetObjResult(interp, Tcl_NewBooleanObj(Odie_IsColinear(x1,y1,x2,y2,x3,y3)));
  return TCL_OK;
}

/*
** This file implements several math and drawing functions used
** to accellerate the IRM gui
*/

static int  odiemath_method_double_to_fuzzy (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double ax;
  if(Tcl_GetDoubleFromObj(interp,objv[1],&ax)) {
    Tcl_SetObjResult(interp,objv[1]);
    return TCL_OK;
  }
  Tcl_SetObjResult(interp,ODIE_NewFuzzyObj(ax));
  return TCL_OK;
}

static int  odiemath_method_fuzzy_compare (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  int result;
  if(objc!=3) {
    Tcl_WrongNumArgs(interp, 1, objv, "avalue bvalue");
  }
  result=ODIE_Fuzzy_Compare_TclObj(objv[1],objv[2]);
  Tcl_SetObjResult(interp,Tcl_NewIntObj(result));
  return TCL_OK;
}

static int  odiemath_method_fuzzy_is_zero (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  int result;
  if(objc!=2) {
    Tcl_WrongNumArgs(interp, 1, objv, "avalue");
  }
  double value;
  if(Tcl_GetDoubleFromObj(NULL,objv[1],&value)) {
    result=0;
  } else {
    result=ODIE_Real_Is_Zero(value);
  }
  Tcl_SetObjResult(interp,Tcl_NewBooleanObj(result));
  return TCL_OK;
}

static int  odiemath_method_grid_hex (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double grid, x, y;
  int gx,gy;
  Tcl_Obj *pResult;
 
  if( objc != 4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "gridsize x y");
    return TCL_ERROR;
  }
  if(Tcl_GetDoubleFromObj(interp,objv[1],&grid)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&x)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&y)) return TCL_ERROR;
  pResult=Tcl_NewObj();
  gy=(int)round(y/grid);
  if(gy%2==1){
    gx=(int)round((x-grid/2)/grid);
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(grid*gx+grid/2));
  } else {
    gx=(int)round(x/grid);
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(grid*gx));
  }
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(grid*gy));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  odiemath_method_grid_square (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double grid;
  double x;
  double y;
  Tcl_Obj *pResult;
  if( objc != 4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "gridsize x y");
    return TCL_ERROR;
  }
  if(Tcl_GetDoubleFromObj(interp,objv[1],&grid)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&x)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&y)) return TCL_ERROR;
  pResult=Tcl_NewObj();
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(grid*round(x/grid)));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(grid*round(y/grid)));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  odiemath_method_line_circle_intersect (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double ax1,ax2,ay1,ay2;
  double bx1,by1,brad;
  double ix,iy;

  if(objc<8) {
    Tcl_WrongNumArgs(interp, 1, objv, "ax1 ay1 ax2 ay2 bx1 by1 brad");
    return TCL_ERROR;
  }

  if( Tcl_GetDoubleFromObj(interp, objv[1], &ax1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[2], &ay1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[3], &ax2) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[4], &ay2) ) return TCL_ERROR;

  if( Tcl_GetDoubleFromObj(interp, objv[5], &bx1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[6], &by1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[7], &brad) ) return TCL_ERROR;
  
  if(ODIE_Math_LineCircleIntersect(ax1,ay1,ax2,ay2,bx1,by1,brad,&ix,&iy)) {
    Tcl_Obj *pResult;
    pResult = Tcl_NewObj();
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ix));
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(iy));
    Tcl_SetObjResult(interp, pResult);
  }
  return TCL_OK;
}

static int  odiemath_method_line_intersect (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double ax1,ax2,ay1,ay2;
  double bx1,bx2,by1,by2;
  double ix,iy;

  if(objc<9) {
    Tcl_WrongNumArgs(interp, 1, objv, "ax1 ay1 ax2 ay2 bx1 by1 bx2 by2");
    return TCL_ERROR;
  }

  if( Tcl_GetDoubleFromObj(interp, objv[1], &ax1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[2], &ay1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[3], &ax2) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[4], &ay2) ) return TCL_ERROR;

  if( Tcl_GetDoubleFromObj(interp, objv[5], &bx1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[6], &by1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[7], &bx2) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[8], &by2) ) return TCL_ERROR;
  
  if(ODIE_Math_LineLineIntersect(ax1,ay1,ax2,ay2,bx1,by1,bx2,by2,&ix,&iy)) {
    Tcl_Obj *pResult;
    pResult = Tcl_NewObj();
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ix));
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(iy));
    Tcl_SetObjResult(interp, pResult);
  }
  return TCL_OK;
}

static int  odiemath_method_line_overlap (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  Tcl_Obj *pResult;
  double ax1,ax2,ay1,ay2;
  double bx1,bx2,by1,by2;

  if(objc<9) {
    Tcl_WrongNumArgs(interp, 1, objv, "ax1 ay1 ax2 ay2 bx1 by1 bx2 by2");
    return TCL_ERROR;
  }

  if( Tcl_GetDoubleFromObj(interp, objv[1], &ax1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[2], &ay1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[3], &ax2) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[4], &ay2) ) return TCL_ERROR;

  if( Tcl_GetDoubleFromObj(interp, objv[5], &bx1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[6], &by1) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[7], &bx2) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[8], &by2) ) return TCL_ERROR;
  
  /*
  ** ignore if the segments connect at endpoints
  if(ax1==bx1 && ay1==by1) return TCL_OK;
  if(ax1==bx2 && ay1==by2) return TCL_OK;
  if(ax2==bx1 && ay2==by1) return TCL_OK;
  if(ax2==bx2 && ay2==by2) return TCL_OK;
  */
    
  pResult = Tcl_NewIntObj(ODIE_Math_LineLineCoincident(ax1,ay1,ax2,ay2,bx1,by1,bx2,by2));
  Tcl_SetObjResult(interp, pResult);

  return TCL_OK;
}

static int  odiemath_method_list_round (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  int i, n;
  double factor;
  if(Tcl_GetDoubleFromObj(interp,objv[1],&factor)) return TCL_ERROR;
  if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;

  Tcl_Obj *pResult=Tcl_NewObj();
  for(i=0;i<n;i++) {
    double thisval;
    Tcl_Obj *pObj;
    Tcl_ListObjIndex(0, objv[2], i, &pObj);
    if(Tcl_GetDoubleFromObj(interp,pObj,&thisval)) {
      Tcl_DecrRefCount(pResult);
      return TCL_ERROR;
    }
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(round(thisval/factor)*factor));
  }
  Tcl_SetObjResult(interp,pResult);
  return TCL_OK;
}

static int  odiemath_method_list_to_int (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  for(i=1;i<objc;i++) {
    double thisval;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&thisval)) return TCL_ERROR;
    Tcl_ListObjAppendElement(interp, pResult, Tcl_NewIntObj((int)thisval));
  }
  Tcl_SetObjResult(interp,pResult);
  return TCL_OK;
}

int Odie_Mathtools_Init(Tcl_Interp *interp) {
  Tcl_Namespace *modPtr;

  modPtr=Tcl_FindNamespace(interp,"odiemath",NULL,TCL_NAMESPACE_ONLY);
  if(!modPtr) {
    modPtr = Tcl_CreateNamespace(interp, "odiemath", NULL, NULL);
  }

  Tcl_CreateObjCommand(interp,"::odiemath::colinear",(Tcl_ObjCmdProc *)odiemath_method_colinear,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::double_to_fuzzy",(Tcl_ObjCmdProc *)odiemath_method_double_to_fuzzy,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::fuzzy_compare",(Tcl_ObjCmdProc *)odiemath_method_fuzzy_compare,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::fuzzy_is_zero",(Tcl_ObjCmdProc *)odiemath_method_fuzzy_is_zero,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::grid_hex",(Tcl_ObjCmdProc *)odiemath_method_grid_hex,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::grid_square",(Tcl_ObjCmdProc *)odiemath_method_grid_square,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::line_circle_intersect",(Tcl_ObjCmdProc *)odiemath_method_line_circle_intersect,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::line_intersect",(Tcl_ObjCmdProc *)odiemath_method_line_intersect,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::line_overlap",(Tcl_ObjCmdProc *)odiemath_method_line_overlap,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::list_round",(Tcl_ObjCmdProc *)odiemath_method_list_round,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::odiemath::list_to_int",(Tcl_ObjCmdProc *)odiemath_method_list_to_int,NULL,NULL);

  Tcl_CreateEnsemble(interp, modPtr->fullName, modPtr, TCL_ENSEMBLE_PREFIX);
  Tcl_Export(interp, modPtr, "[a-z]*", 1);
  return TCL_OK;
}

/*
** This widget translates 3-D coordinates onto a flat canvas by splitting
** the 3-D space into layers and stacking the layers on the canvas.
**
** The layers are decks of the ship.  The highest layer (or deck) is drawn
** at the top of the page.  The next layer down is drawn below the top layer.
** and so forth down the canvas.  In other words, the 3D object is drawn
** by showing a set of 2D slices where each slice is viewed from above.
**
** The original 3D coordinates are called "actual" coordinates.  When
** translated into the 2D canvas they are called "canvas" coordinates.
**
** The actual coordinate system is right-handed.  The X axis increases to
** the right.  The Y axis increases going up.  The Z axis comes out of the
** page at the viewer.  The canvas coordinate system is left-handed.  The
** X axis increase to the right but the Y axis increases going down.
**
** A plotter is a object with methods.  The details of the available
** methods and what each does are described in comments before the
** implementation of each method.
*/
#include "odieInt.h"
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h>

/*
** A plotter is an instance of the following structure.
*/
typedef struct Plotter Plotter;
struct Plotter {
  double rZoom;      /* Multiply canvas coord by this to get actual coord */
  double rXOffset;   /* X-Shift amount */
  double rYOffset;   /* Y-Shift amount */
};

/*
** This routine is called when a plotter is deleted.  All the memory and
** other resources allocated by this plotter is recovered.
*/
static void destroyPlotter(void *pArg){
  Plotter *p = (Plotter*)pArg;
  Odie_Free((char *)p);
}

static inline double xCanvasToActual(Plotter *p,double cx){
  return (cx+p->rXOffset)*p->rZoom;
}

static inline double yCanvasToActual(Plotter *p,double cy){
  return -1.0*(cy+p->rYOffset)*p->rZoom;
}

/*
** Convert a Y coordinate from actual to canvas coordinates for a
** given deck.
*/
static inline double xActualToCanvas(Plotter *p,double ax){
  return (ax/p->rZoom)-p->rXOffset;
}

/*
** Convert a Y coordinate from actual to canvas coordinates for a
** given deck.
*/
static inline double yActualToCanvas(Plotter *p,double ay){
  return -1.0*(ay/p->rZoom)-p->rYOffset;
}

/*
** This routine runs when a method is executed against a plotter
*/
static int plotterMethodProc(
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Plotter *p = (Plotter*)pArg;

#if 0
  /* For debugging....
  ** Print each wallset command before it is executed.
  */
  { int i;
    for(i=0; i<objc; i++){
      printf("%s%c", Tcl_GetStringFromObj(objv[i], 0), i<objc-1 ? ' ' : '\n');
    }
  }
#endif

  /* The following bit of magic implements a switch() statement that
  ** invokes the correct code based on the value of objv[1].  The
  ** mktclopts.tcl script scans this source file looking for "case"
  ** statements then generates the "wallset.h" include file that contains
  ** all of the code necessary to implement the switch.
  **
  ** In this way, we can add new methods to the command simply by adding
  ** new cases within the switch.  All the related switch code is
  ** regenerated automatically.
  */
  {
#include "plotter_cases.h"
  {
  /*
  ** tclmethod:  PLOTTER actualcoords CANVAS-COORD-LIST
  ** title:   Convert from canvas to actual coordinate space
  */
  case PLOTTER_ACTUALCOORDS: {
    int i, n;
    Tcl_Obj *pResult;
    MATOBJ *M;
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "CANVAS-COORD-LIST");
      return TCL_ERROR;
    }
    if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;
    if( n%2!=0 ){
      Tcl_AppendResult(interp, "coordinate list must contain a multiple "
         "of 2 values", 0);
      return TCL_ERROR;
    }

    pResult = Tcl_NewObj();
    for(i=0; i<n-1; i+=2){
      double cx, cy;
      Tcl_Obj *pObj;

      Tcl_ListObjIndex(0, objv[2], i, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &cx) ) break;
      Tcl_ListObjIndex(0, objv[2], i+1, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &cy) ) break;
      M=Matrix_NewObj();
      Matrix_Alloc(M,MATFORM_vectorxy);

      /* Original Formula 
      ** ax = cx*p->rZoom - pS->rXShift;
      ** ay = pS->mxY + (pS->mnY-pS->mxY)*(cy-pS->top)/(pS->btm-pS->top);
      */
      *(M->matrix+X_IDX)=xCanvasToActual(p,cx);
      *(M->matrix+Y_IDX)=yCanvasToActual(p,cy);
      
      Tcl_ListObjAppendElement(interp, pResult, Matrix_To_TclObj(M));
    }
    if( i<n-1 ){
      Tcl_DecrRefCount(pResult);
      return TCL_ERROR;
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }
 
  /*
  ** tclmethod:  PLOTTER canvascoords VECTOR-LIST ?xvar? ?yvar?
  ** title:   Convert from actual (a list of vectors) to canvas coordinate space
  **
  */
  case PLOTTER_CANVASCOORDS: {
    int i,n,m;
    Tcl_Obj *pResult=NULL;
    double cx, cy, ax, ay;
    Tcl_Obj *pObj;
    int error=0;
    int singlearg=0;

    if( objc!=3 && objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "VECTOR-LIST ?xvar yvar?");
      return TCL_ERROR;
    }
    
    
    if(objv[2]->typePtr==&matrix_tclobjtype) {
      singlearg=1;
    } else {
      if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;
      Tcl_ListObjIndex(0, objv[2], 0, &pObj);
      if( Tcl_ListObjLength(interp, pObj, &m) ) return TCL_ERROR;
      if(m==1) {
        singlearg=1;
      }
    }
    
    if(singlearg) {
      /*
      ** Accept a single vector as an argument
      ** Do this to ensure we don't interpret the
      ** value as a list
      */
      MATOBJ *M;
      M=Odie_GetMatrixFromTclObj(interp,objv[2],MATFORM_vectorxy);
      if(!M) return TCL_ERROR;
      ax=*(M->matrix+X_IDX);
      ay=*(M->matrix+Y_IDX);
      cx = xActualToCanvas(p,ax);       
      cy = yActualToCanvas(p,ay);
      if(objc==5) {
        Tcl_ObjSetVar2(interp,objv[3],NULL,Tcl_NewDoubleObj(cx),0);
        Tcl_ObjSetVar2(interp,objv[4],NULL,Tcl_NewDoubleObj(cy),0);
      }
      pResult = Tcl_NewObj();
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cx));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cy));
      Tcl_SetObjResult(interp, pResult);  
      return TCL_OK;
    }

    if(objc != 5) {
      pResult = Tcl_NewObj();
    }

    for(i=0; i<n; i++){
      MATOBJ *M;
      Tcl_ListObjIndex(0, objv[2], i, &pObj);
      M=Odie_GetMatrixFromTclObj(interp,pObj,MATFORM_vectorxy);
      if(!M) return TCL_ERROR;
      ax=*(M->matrix+X_IDX);
      ay=*(M->matrix+Y_IDX);
      cx = xActualToCanvas(p,ax);       
      cy = yActualToCanvas(p,ay);
      if(objc==5) {
        Tcl_ObjSetVar2(interp,objv[3],NULL,Tcl_NewDoubleObj(cx),0);
        Tcl_ObjSetVar2(interp,objv[4],NULL,Tcl_NewDoubleObj(cy),0);
        return TCL_OK;
      } else {    
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cx));
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cy));
      }
    }
    if(error) {
      if(pResult) {
        Tcl_DecrRefCount(pResult);
      }
      return TCL_ERROR;
    }
    if( i<n-3 ){
        Tcl_AppendResult(interp, "Did not reach the end", 0);
      Tcl_DecrRefCount(pResult);
      return TCL_ERROR;
    }

    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  PLOTTER destroy
  ** title:   Destroy this plotter
  */
  case PLOTTER_DESTROY: {
    Tcl_DeleteCommand(interp,Tcl_GetString(objv[0]));
    break;
  }

  /*
  ** tclmethod:  PLOTTER objinfo obj
  ** title:   Return TK to draw a line on a canvas
  */
  case PLOTTER_OBJINFO: {
    Tcl_Obj *tmp;
    tmp = objv[2];
    printf("INFO: %s Ref: %d Type: %p \n", Tcl_GetStringFromObj(tmp, NULL),
        tmp->refCount, tmp->typePtr);
    fflush (stdout);
    break;
  }
  
  /*
  ** tclmethod:  PLOTTER centerset zoom width height
  ** title:   Set all settings for plotter in one go
  ** description: Sets the center of the screen based on the width
  ** and height (0,0 = width/2 height/2)
  */
  case PLOTTER_CENTERSET: {
    double rZoom,rXOffset,rYOffset;

    if(objc!=5) {
      printf("%d\n",objc);
      Tcl_WrongNumArgs(interp, 2, objv, "ZOOM XOFFSET YOFFSET");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &rZoom) ) return TCL_ERROR;
    p->rZoom = rZoom;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &rXOffset) ) return TCL_ERROR;
    p->rXOffset = -rXOffset/2.0;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &rYOffset) ) return TCL_ERROR;
    p->rYOffset = -rYOffset/2.0;
    return TCL_OK;
  }
  
  /*
  ** tclmethod:  PLOTTER xoffset ?AMT?
  ** title:   Change the X-Offset
  */
  case PLOTTER_XOFFSET: {
    double rXOffset;
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==2 ){
      Tcl_SetObjResult(interp, Tcl_NewDoubleObj(p->rXOffset));
      return TCL_OK;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &rXOffset) ) return TCL_ERROR;
    p->rXOffset = rXOffset;
    break;
  }

  /*
  ** tclmethod:  PLOTTER yoffset ?AMT?
  ** title:   Change the Y-Offset
  */
  case PLOTTER_YOFFSET: {
    double rYOffset;
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==2 ){
      Tcl_SetObjResult(interp, Tcl_NewDoubleObj(p->rYOffset));
      return TCL_OK;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &rYOffset) ) return TCL_ERROR;
    p->rYOffset = rYOffset;
    break;
  }

  /*
  ** tclmethod:  PLOTTER zoom ?ZOOM?
  ** title:   Query or change the zoom factor.
  */
  case PLOTTER_ZOOM: {
    Tcl_Obj *pResult;
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==3 ){
      double r;
      if( Tcl_GetDoubleFromObj(interp, objv[2], &r) ) return TCL_ERROR;
      p->rZoom = r;
    }
    pResult = Tcl_NewDoubleObj(p->rZoom);
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /* End of the command methods.  The brackets that follow terminate the
  ** automatically generated switch.
  ****************************************************************************/
  }
  }
  return TCL_OK;
}

/*
** tclcmd: plotter PLOTTER
** title: creates a plotter object
** This routine runs when the "plotter" command is invoked to create a
** new plotter.
*/
int Odie_PlotterCreateProc(
  void *NotUsed,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  char *zCmd;
  Plotter *p;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "PLOTTER");
    return TCL_ERROR;
  }
  zCmd = Tcl_GetStringFromObj(objv[1], 0);
  p = (Plotter *)Odie_Alloc( sizeof(*p) );
  p->rZoom = 1.0;
  p->rXOffset = 0.0;
  p->rYOffset = 0.0;
  Tcl_CreateObjCommand(interp, zCmd, plotterMethodProc, p, destroyPlotter);
  return TCL_OK;
}

/*** Automatically Generated Header File - Do Not Edit ***/
  const static char *PLOTTER_strs[] = {
    "actualcoords",       "canvascoords",      "centerset",
    "destroy",            "objinfo",           "xoffset",
    "yoffset",            "zoom",              0
  };
  enum PLOTTER_enum {
    PLOTTER_ACTUALCOORDS, PLOTTER_CANVASCOORDS,PLOTTER_CENTERSET,
    PLOTTER_DESTROY,      PLOTTER_OBJINFO,     PLOTTER_XOFFSET,
    PLOTTER_YOFFSET,      PLOTTER_ZOOM,        
  };
 int index;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "METHOD ?ARG ...?");
    return TCL_ERROR;
  }
  if( Tcl_GetIndexFromObj(interp, objv[1], PLOTTER_strs, "option", 0, &index)){
    return TCL_ERROR;
  }
  switch( (enum PLOTTER_enum)index )

/*
** This file is machine generated. Changes will
** be overwritten on the next run of cstruct.tcl
*/
#include "odieInt.h"

/*
** This file implements a TCL object used for tracking polygons.  A
** single new TCL command named "poly" is defined.  This command
** has methods for creating, deleting, and taking the intersection
** of 2-D polygons.  There are comments on the implementation of each
** method to describe what the method does.
**
** This module was originally developed to aid in computing the
** shared surface area between two compartments on separate decks.
** The shared surface area is needed in initializing the fire model
** since heat conduction between the two compartments is proportional
** to their shared area.
*/
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <math.h>

/*
** Return the distance between two points
*/
static inline double dist(double x0, double y0, double x1, double y1){
  double dx = x1 - x0;
  double dy = y1 - y0;
  return sqrt(dx*dx + dy*dy);
}

/*
** Return -1, 0, or 1 if the point x,y is outside, on, or within
** the polygon p.
*/
static inline int within(Poly *p, double x, double y){
  int res, i;
  res = -1;
  for(i=0; i<p->nVertex-1; i++){
    double x0, y0, x1, y1, yP;
    x0 = p->v[i][X_IDX];
    y0 = p->v[i][Y_IDX];
    x1 = p->v[i+1][X_IDX];
    y1 = p->v[i+1][Y_IDX];
    if( x0==x1 ){
      if( x0==x && ((y0<=y && y1>=y) || (y1<=y && y0>=y)) ){
        res = 0;
        break;
      }
      continue;
    }
    if( x0>x1 ){
      int t = x0;
      x0 = x1;
      x1 = t;
      t = y0;
      y0 = y1;
      y1 = t;
    }
    if( x>=x1 || x<x0 ) continue;
    yP = y1 - (x1-x)*(y1-y0)/(x1-x0);
    if( yP == y ){ res = 0; break; }
    if( yP > y ){ res = -res; }
  }
  return res;
}

int Odie_PolygonComputeArea(Tcl_Interp *interp,Poly *p) {
  double area=0.0;
  double areacomp=0.0;
  int i;

  if( p->v[p->nVertex-1][X_IDX]!=p->v[0][X_IDX] || p->v[p->nVertex-1][Y_IDX]!=p->v[0][Y_IDX] ){
    p->v[p->nVertex][X_IDX] = p->v[0][X_IDX];
    p->v[p->nVertex][Y_IDX] = p->v[0][Y_IDX];
    p->nVertex++;
  }

  for(i=0; i<p->nVertex-1; i++){
    area += 0.5*(p->v[i][Y_IDX] + p->v[i+1][Y_IDX])*(p->v[i+1][X_IDX] - p->v[i][X_IDX]);
  }
  if( area<0.0 ){
    int b, e;
    for(b=0, e=p->nVertex-1; b<e; b++, e--){
      double t;
      t = p->v[b][X_IDX];
      p->v[b][X_IDX] = p->v[e][X_IDX];
      p->v[e][X_IDX] = t;
      t = p->v[b][Y_IDX];
      p->v[b][Y_IDX] = p->v[e][Y_IDX];
      p->v[e][Y_IDX] = t;
    }
    area = -area;
  }
  p->area = area;
  p->bbox.l = p->bbox.r = p->v[0][X_IDX];
  p->bbox.t = p->bbox.b = p->v[0][Y_IDX];
  for(i=1; i<p->nVertex-1; i++){
    double x, y;
    x = p->v[i][X_IDX];
    if( x<p->bbox.l ) p->bbox.l = x;
    if( x>p->bbox.r ) p->bbox.r = x;
    y = p->v[i][Y_IDX];
    if( y>p->bbox.t ) p->bbox.t = y;
    if( y<p->bbox.b ) p->bbox.b = y;
  }
  areacomp=(p->bbox.r - p->bbox.l)*(p->bbox.t-p->bbox.b)*1.00001;
  
  if(area<=areacomp) {
    return TCL_OK;
  } else {
    char errstr[256];
    sprintf(errstr,"Area: %g Calculated: %g\n",area,areacomp);
    Tcl_AppendResult(interp, "Area of polygon wonky ", errstr, 0);
    return TCL_ERROR;
  }
}

static int  polygon_method_create (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N");
    return TCL_ERROR;
  }
  int isnew;
  Poly *p;
  if( Odie_GetPolygonFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  if(isnew) {
    Tcl_SetObjResult(interp, Odie_NewPolygonObj(p));
  } else {
    Tcl_SetObjResult(interp, objv[1]);
  }
  return TCL_OK;
}

static int  polygon_method_simplify (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N");
    return TCL_ERROR;
  }
  Poly *pPoly,*pNewPoly;
  int i,isnew;
  int colinear;
  double ax,ay,bx,by,cx,cy;
  if( Odie_GetPolygonFromObj(interp, objv[1], &pPoly, &isnew) ) return TCL_ERROR;

  pNewPoly=(Poly *)Odie_Alloc(sizeof(*pNewPoly)+(pPoly->nVertex+2)*sizeof(pNewPoly->v[0]));
  pNewPoly->nVertex=0;
  
  ax=pPoly->v[pPoly->nVertex-1][X_IDX];
  ay=pPoly->v[pPoly->nVertex-1][Y_IDX];
  bx=pPoly->v[0][X_IDX];
  by=pPoly->v[0][Y_IDX];
  if(ax==bx && ay==by) {
    ax=pPoly->v[pPoly->nVertex-2][X_IDX];
    ay=pPoly->v[pPoly->nVertex-2][Y_IDX];
  }
  for(i=1;i<pPoly->nVertex;i++) {
    cx=pPoly->v[i][X_IDX];
    cy=pPoly->v[i][Y_IDX];
    colinear=Odie_IsColinear(ax,ay,bx,by,cx,cy);
    if(!colinear) {
      pNewPoly->v[pNewPoly->nVertex][X_IDX]=bx;
      pNewPoly->v[pNewPoly->nVertex][Y_IDX]=by;
      pNewPoly->nVertex++;
    }
    ax=bx;
    ay=by;
    bx=cx;
    by=cy;
  }
  cx=pPoly->v[0][X_IDX];
  cy=pPoly->v[0][Y_IDX];
  colinear=Odie_IsColinear(ax,ay,bx,by,cx,cy);
  if(!Odie_IsColinear(ax,ay,bx,by,cx,cy)) {
    pNewPoly->v[pNewPoly->nVertex][X_IDX]=bx;
    pNewPoly->v[pNewPoly->nVertex][Y_IDX]=by;
    pNewPoly->nVertex++;
  }
  if( pNewPoly->v[pNewPoly->nVertex-1][X_IDX]!=pNewPoly->v[0][X_IDX] || pNewPoly->v[pNewPoly->nVertex-1][Y_IDX]!=pNewPoly->v[0][Y_IDX] ){
    pNewPoly->v[pNewPoly->nVertex][X_IDX] = pNewPoly->v[0][X_IDX];
    pNewPoly->v[pNewPoly->nVertex][Y_IDX] = pNewPoly->v[0][Y_IDX];
    pNewPoly->nVertex++;
  }
  Odie_PolygonComputeArea(interp,pNewPoly);
  Tcl_SetObjResult(interp, Odie_NewPolygonObj(pNewPoly));
  if(isnew) Odie_Free((char *)pPoly);
  return TCL_OK;
}

static int  polygon_method_area (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  Poly *p;
  int isnew;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N ?N...?");
    return TCL_ERROR;
  }
  if( Odie_GetPolygonFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  double area=p->area;
  if(isnew) Odie_Free((char *)p);
  int i;
  for(i=2;i<objc;i++) {
    if( Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew) ) return TCL_ERROR;
    area+=p->area;
    if(isnew) Odie_Free((char *)p);
  }
  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(area));
  return TCL_OK;
}

static int  polygon_method_bbox (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*    poly bbox N
  **
  ** Return the bounding box for a polygon
  */
  Poly *p;
  int isnew;
  Tcl_Obj *pResult;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N ?N...?");
    return TCL_ERROR;
  }
  if( Odie_GetPolygonFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  double left=p->bbox.l;
  double top=p->bbox.t;
  double right=p->bbox.r;
  double bottom=p->bbox.b;
  int i;
  if(isnew) Odie_Free((char *)p);
  for(i=2;i<objc;i++) {
    if( Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew) ) return TCL_ERROR;
    if(p->bbox.l < left) left=p->bbox.l;
    if(p->bbox.t > top)  top=p->bbox.t;
    if(p->bbox.r > right) right=p->bbox.r;
    if(p->bbox.b < bottom)  bottom=p->bbox.b;
    if(isnew) Odie_Free((char *)p);
  }
  pResult = Tcl_NewObj();
  Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(left));
  Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(top));
  Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(right));
  Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(bottom));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_info (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*    poly info N
  **
  ** Return the coordinates for a polygon
  */
  Tcl_Obj *pResult;
  Poly *p;
  int i,isnew;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N");
    return TCL_ERROR;
  }
  if( Odie_GetPolygonFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  pResult = Tcl_NewObj();
  for(i=0; i<p->nVertex-1; i++){
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(p->v[i][X_IDX]));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(p->v[i][Y_IDX]));
  }
  if(isnew) Odie_Free((char *)p);
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_intersect (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*    poly intersect N1 N2
  **
  ** Return a list of 3 elements where the first element is the
  ** area of intersection between polygons N1 and N2 and the remaining
  ** 3 elements are the X and Y coordinates of a point within both
  ** polygons. 
  **
  ** The current implementation returns an approximation.  We might
  ** change it to compute the exact intersection later.
  */    
  Poly *p1, *p2;
  int isnew1,isnew2;
  double area;
  double xInside = 0.0, yInside = 0.0;
  Tcl_Obj *pResult;
  if( objc!=3 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N1 N2");
    return TCL_ERROR;
  }
  if( Odie_GetPolygonFromObj(interp, objv[1], &p1, &isnew1) ) return TCL_ERROR;
  if( Odie_GetPolygonFromObj(interp, objv[2], &p2, &isnew2) ) {
    if(isnew1) Odie_Free((char *)p1);
    return TCL_ERROR;
  }
  if( p1->bbox.r<=p2->bbox.l || p1->bbox.l>=p2->bbox.r
                || p1->bbox.t<=p2->bbox.b  || p1->bbox.b>=p2->bbox.t ){
    area = 0.0;
  }else if( p1->area==0.0 || p2->area==0.0 ){
    area = 0.0;
  }else{
    double x0, y0, x1, y1, dx, dy, xP, yP, xC, yC;
    int i, j, cnt;
    int score, bestScore;
    static const int n = 50;
    char hit[50][50];

    /* Compute the overlap of the bounding boxes of the two polygons. */
    x0 = p1->bbox.l < p2->bbox.l ? p2->bbox.l : p1->bbox.l;
    y0 = p1->bbox.t > p2->bbox.t ? p2->bbox.t : p1->bbox.t;
    x1 = p1->bbox.r > p2->bbox.r ? p2->bbox.r : p1->bbox.r;
    y1 = p1->bbox.b < p2->bbox.b ? p2->bbox.b : p1->bbox.b;

    /* Divide the intersection of the bounding boxes into a n-by-n grid
    ** and count the number of elements in this grid whose centers fall
    ** within both polygons.  This will be our approximation for the
    ** intersection of the polygons themselves.
    */
    dx = (x1-x0)/n;
    dy = (y1-y0)/n;
    cnt = 0;
    xC = yC = 0.0;
    for(i=0; i<n; i++){
      xP = x0 + dx*(i+0.5);
      for(j=0; j<n; j++){
        yP = y0 + dy*(j+0.5);
        if( within(p1, xP, yP)>0 && within(p2, xP, yP)>0 ){
          cnt++;
          hit[i][j] = 1;
          xC += xP;
          yC += yP;
        }else{
          hit[i][j] = 0;
        }
      }
    }

    /* We need to find a good approximation for the center of the
    ** overlap.  Begin by computing the center of mass for the
    ** overlapping region.  Then find the point inside the intersection
    ** that is nearest the center of mass.
    */
    if( cnt>0 ){
      area = cnt*(x1-x0)*(y0-y1)/(n*n);
      xC /= cnt;
      yC /= cnt;
      bestScore = -1.0;
      for(i=0; i<n; i++){
        xP = x0 + dx*(i+0.5);
        for(j=0; j<n; j++){
          if( !hit[i][j] ) continue;
          yP = y0 + dy*(j+0.5);
          score = dist(xP,yP,xC,yC);
          if( score<bestScore || bestScore<0.0 ){
            xInside = xP;
            yInside = yP;
            bestScore = score;
          }
        }
      }
    } else {
      area=0.0;
    }
  }
  if(isnew1) Odie_Free((char *)p1);
  if(isnew2) Odie_Free((char *)p2);

  pResult = Tcl_NewObj();
  Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(area));
  Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(xInside));
  Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(yInside));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_within (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*    poly within X Y N ?N...?
  **
  ** Return -1, 0, or +1 if the point X,Y is outside, on, or inside
  ** polygon N.
  */
  Poly *p;
  int res,isnew;
  double x, y;
  if( objc<4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "X Y ?N...?");
    return TCL_ERROR;
  }
  if( Tcl_GetDoubleFromObj(interp, objv[1], &x) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, objv[2], &y) ) return TCL_ERROR;
  if(objc==4) {
    if( Odie_GetPolygonFromObj(interp, objv[3], &p, &isnew) ) return TCL_ERROR;
    res = within(p, x, y);
    if(isnew) Odie_Free((char *)p);
    Tcl_SetObjResult(interp, Tcl_NewIntObj(res));
    return TCL_OK;      
  } else {
    int i;
    for(i=3;i<objc;i++) {
      if( Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew) ) return TCL_ERROR;
      res = within(p, x, y);
      if(isnew) Odie_Free((char *)p);
      if(res>=0) {
        Tcl_SetObjResult(interp, Tcl_NewBooleanObj(1));
        return TCL_OK;      
      }
    }
    Tcl_SetObjResult(interp, Tcl_NewBooleanObj(0));
    return TCL_OK;
  }
}

static int  polygon_method_segments (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  Tcl_Obj *pResult;
  int isnew;
  Poly *p;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N");
    return TCL_ERROR;
  }
  if( Odie_GetPolygonFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  int i;
  double px,py;
  pResult=Tcl_NewObj();
  px=p->v[0][X_IDX];
  py=p->v[0][Y_IDX];
  for(i=1; i<p->nVertex-1; i++){
    Tcl_Obj *segment=Tcl_NewObj();
    Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(px));
    Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(py));
    Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(p->v[i][X_IDX]));
    Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(p->v[i][Y_IDX]));
    Tcl_ListObjAppendElement(interp,pResult, segment);
    px=p->v[i][X_IDX];
    py=p->v[i][Y_IDX];
  }
  Tcl_Obj *segment=Tcl_NewObj();
  Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(px));
  Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(py));
  Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(p->v[0][X_IDX]));
  Tcl_ListObjAppendElement(interp,segment, Tcl_NewDoubleObj(p->v[0][Y_IDX]));
  Tcl_ListObjAppendElement(interp,pResult, segment);
  Tcl_SetObjResult(interp, pResult);
  if(isnew) Odie_Free((char *)p);
  return TCL_OK;
}

static int  polygon_method_rectangle (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double cx, cy, radx,rady;
  Tcl_Obj *pResult=Tcl_NewObj();

  if( objc != 5 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  radx=radx/2.0;
  rady=rady/2.0;
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}


static int  polygon_method_vector_place (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*
  ** Apply Matrices
  */
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double zoom;
  double matA[6] = {1.0,0.0,0.0,1.0,0.0,0.0};
  double centerx,centery,normalx,normaly,angle;

  if( objc < 8 ){
      Tcl_WrongNumArgs(interp, 1, objv, "zoom centerx centery normalx normaly x1 y1 ?x2 y2?...");
      return TCL_ERROR;
  }   
  if(Tcl_GetDoubleFromObj(interp,objv[1],&zoom)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&centerx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&centery)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&normalx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[5],&normaly)) return TCL_ERROR;

  angle=atan2(normaly,normalx);
  matA[0]=cos(angle);
  matA[1]=sin(angle);
  matA[2]=-sin(angle);
  matA[3]=cos(angle);
  matA[4]=0.0;
  matA[5]=0.0;

  
  for(i=6;i<objc;i+=2) {
      double x,y,sx,sy,newx,newy;
      if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
      if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
      
      sx=(x/zoom);
      sy=(y/zoom);
      newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
      newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newx));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newy));
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}


static int  polygon_method_hexagon (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  int i,flip=0;
  double cx, cy, radx,rady;

  Tcl_Obj *pResult=Tcl_NewObj();
  double coords[7][2]= {
    {1.00, 0.00} , {0.50, M_SQRT3_2} , 
    {-0.50, M_SQRT3_2} , {-1.00, -0.00} , 
    {-0.50, -M_SQRT3_2} , {0.50, -M_SQRT3_2},
    {1.00, 0.00} 
  };
  if( objc != 5 && objc != 6){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy ?flip?");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  if(objc==6) {
    if(Tcl_GetBooleanFromObj(interp,objv[5],&flip)) return TCL_ERROR;
  }
  radx=radx/2.0;
  rady=rady/2.0;

  for(i=0;i<6;i++) {
    if(flip) {
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i+1][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i+1][0]));
    } else {
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i+1][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i+1][1]));
    }
  }

  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_poly_place (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*
  ** Apply Matrices
  */
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double zoom;
  double matA[6] = {1.0,0.0,0.0,1.0,0.0,0.0};
  double centerx,centery,normalx,normaly,angle;

  if( objc < 8 ){
      Tcl_WrongNumArgs(interp, 1, objv, "zoom centerx centery normalx normaly x1 y1 ?x2 y2?...");
      return TCL_ERROR;
  }   
  if(Tcl_GetDoubleFromObj(interp,objv[1],&zoom)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&centerx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&centery)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&normalx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[5],&normaly)) return TCL_ERROR;

  angle=atan2(normaly,normalx);
  matA[0]=cos(angle);
  matA[1]=sin(angle);
  matA[2]=-sin(angle);
  matA[3]=cos(angle);
  matA[4]=0.0;
  matA[5]=0.0;
  double startx,starty,prevx,prevy;

  i=6;
  {
    double x,y,sx,sy,newx,newy;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
    
    sx=(x/zoom);
    sy=(y/zoom);
    newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
    newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

    startx=newx;
    starty=newy;
    prevx=newx;
    prevy=newy;
  }
  
  for(i=8;i<objc;i+=2) {
    double x,y,sx,sy,newx,newy;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
    
    sx=(x/zoom);
    sy=(y/zoom);
    newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
    newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevy));  
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newy));       

    prevx=newx;
    prevy=newy;
  }
  if(startx != prevx && starty!= prevy) {
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevy));  
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(startx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(starty));  
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_drawobj_orientation (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  Tcl_Obj *temp;
  int len;
  double nx=100,ny=0;
  if( objc !=4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "orientation nxvar nyvar");
    return TCL_ERROR;
  }
  
  if(Tcl_ListObjLength(interp,objv[1],&len)) return TCL_ERROR;
  if(len>0) {
    if(Tcl_ListObjIndex(interp, objv[1], 0, &temp)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,temp,&nx)) return TCL_ERROR;
  }
  if(len>1) {
    if(Tcl_ListObjIndex(interp, objv[1], 1, &temp)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,temp,&ny)) return TCL_ERROR;
  }
  Tcl_ObjSetVar2(interp,objv[2],NULL,Tcl_NewDoubleObj(nx),0);
  Tcl_ObjSetVar2(interp,objv[3],NULL,Tcl_NewDoubleObj(ny),0);
  return TCL_OK;
}

static int  polygon_method_corners (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double cx, cy, radx,rady;
  
  if( objc != 5 && objc != 9 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy ?x0var y0var x1var y1var?");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  if (objc == 5) {
    Tcl_Obj *pResult=Tcl_NewObj();
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));
    Tcl_SetObjResult(interp, pResult);
    return TCL_OK;
  }
  /*
    Replaces
    set x0 [expr {$cx-$d}]
    set y0 [expr {$cy-$d}]
    set x1 [expr {$cx+$d}]
    set y1 [expr {$cy+$d}]
  */
  
  Tcl_ObjSetVar2(interp,objv[5],NULL,Tcl_NewDoubleObj(cx+radx),0);
  Tcl_ObjSetVar2(interp,objv[6],NULL,Tcl_NewDoubleObj(cy-rady),0);
  Tcl_ObjSetVar2(interp,objv[7],NULL,Tcl_NewDoubleObj(cx-radx),0);
  Tcl_ObjSetVar2(interp,objv[8],NULL,Tcl_NewDoubleObj(cy+rady),0);
  
  return TCL_OK; 
}



static int  polygon_method_hexgrid_create (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {

  /*    poly hexgrid
  **
  ** Reduce the polygons to a series
  ** of grid coordinates
  */
  Tcl_Obj *pResult;
  Poly *p;
  int isnew;
  pResult = Tcl_NewObj();
  double gridsize=500.0;
  double x,y;

  if(objc < 3) {
    Tcl_WrongNumArgs(interp, 1, objv, "GRIDSIZE ?POLY ...?");
    return TCL_ERROR;
  }
  if( Tcl_GetDoubleFromObj(interp, objv[1], &gridsize) ) return TCL_ERROR;
  
  if( Odie_GetPolygonFromObj(interp, objv[2], &p, &isnew) ) return TCL_ERROR;
  double left=p->bbox.l;
  double top=p->bbox.t;
  double right=p->bbox.r;
  double bottom=p->bbox.b;
  int i;
  if(isnew) Odie_Free((char *)p);
  for(i=3;i<objc;i++) {
    if( Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew) ) return TCL_ERROR;
    if(p->bbox.l < left) left=p->bbox.l;
    if(p->bbox.t > top)  top=p->bbox.t;
    if(p->bbox.r > right) right=p->bbox.r;
    if(p->bbox.b < bottom)  bottom=p->bbox.b;
    if(isnew) Odie_Free((char *)p);
  }
  
  pResult = Tcl_NewObj();
  left-=gridsize;
  top+=gridsize;
  right+=gridsize;
  bottom-=gridsize;
  int row=0;
  for(y=bottom;y<=top;y+=gridsize) {
    double lstartx=left;
    double gy=floor(y/gridsize)*gridsize;
    row++;
    if(row%2==1) {
      lstartx-=gridsize/2;
    }
    for(x=lstartx;x<=right;x+=gridsize) {
      double gx=floor(x/gridsize)*gridsize;
      int found=0;
      for(i=2;i<objc;i++) {
        int isnew,isWithin;
        Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew);
        isWithin=within(p,gx,gy);
        if(isnew) Odie_Free((char *)p);
        if(isWithin>0) {
          found=1;
          break;
        }
      }
      if(found) {
        Tcl_Obj *coord=Tcl_NewObj();
        Tcl_ListObjAppendElement(interp, coord, Tcl_NewDoubleObj(gx));
        Tcl_ListObjAppendElement(interp, coord, Tcl_NewDoubleObj(gy));
        Tcl_ListObjAppendElement(interp, pResult, coord);
      }
    }
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}


static int  polygon_method_squaregrid_create (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {

  /*    poly hexgrid
  **
  ** Reduce the polygons to a series
  ** of grid coordinates
  */
  Tcl_Obj *pResult;
  Poly *p;
  int isnew;
  pResult = Tcl_NewObj();
  double gridsize=500.0;
  double x,y;
  if(objc < 3) {
    Tcl_WrongNumArgs(interp, 1, objv, "GRIDSIZE ?POLY ...?");
    return TCL_ERROR;
  }
  if( Tcl_GetDoubleFromObj(interp, objv[1], &gridsize) ) return TCL_ERROR;
  
  if( Odie_GetPolygonFromObj(interp, objv[2], &p, &isnew) ) return TCL_ERROR;
  double left=p->bbox.l;
  double top=p->bbox.t;
  double right=p->bbox.r;
  double bottom=p->bbox.b;
  if(isnew) Odie_Free((char *)p);

  int i;
  for(i=3;i<objc;i++) {
    if( Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew) ) return TCL_ERROR;
    if(p->bbox.l < left) left=p->bbox.l;
    if(p->bbox.t > top)  top=p->bbox.t;
    if(p->bbox.r > right) right=p->bbox.r;
    if(p->bbox.b < bottom)  bottom=p->bbox.b;
    if(isnew) Odie_Free((char *)p);
  }
  
  pResult = Tcl_NewObj();
  left-=gridsize;
  top+=gridsize;
  right+=gridsize;
  bottom-=gridsize;

  for(y=bottom;y<=top;y+=gridsize) {
    double gy=floor(y/gridsize)*gridsize;
    for(x=left;x<=right;x+=gridsize) {
      double gx=floor(x/gridsize)*gridsize;
      int found=0;
      for(i=2;i<objc;i++) {
        int isnew,isWithin;
        Odie_GetPolygonFromObj(interp, objv[i], &p, &isnew);
        isWithin=within(p,gx,gy);
        if(isnew) Odie_Free((char *)p);
        if(isWithin>0) {
          found=1;
          break;
        }
      }
      if(found) {
        Tcl_Obj *coord=Tcl_NewObj();
        Tcl_ListObjAppendElement(interp, coord, Tcl_NewDoubleObj(gx));
        Tcl_ListObjAppendElement(interp, coord, Tcl_NewDoubleObj(gy));
        Tcl_ListObjAppendElement(interp, pResult, coord);
      }
    }
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  polygon_method_grid_nearest (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  return TCL_OK;
}

int Odie_Polygon_Init(Tcl_Interp *interp) {
  Tcl_Namespace *modPtr;
  
  modPtr=Tcl_FindNamespace(interp,"polygon",NULL,TCL_NAMESPACE_ONLY);
  if(!modPtr) {
    modPtr = Tcl_CreateNamespace(interp, "polygon", NULL, NULL);
  }
  
  Tcl_CreateObjCommand(interp,"::polygon::area",(Tcl_ObjCmdProc *)polygon_method_area,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::bbox",(Tcl_ObjCmdProc *)polygon_method_bbox,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::corners",(Tcl_ObjCmdProc *)polygon_method_corners,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::create",(Tcl_ObjCmdProc *)polygon_method_create,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::drawobj_orientation",(Tcl_ObjCmdProc *)polygon_method_drawobj_orientation,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::grid_nearest",(Tcl_ObjCmdProc *)polygon_method_grid_nearest,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::hexgrid_create",(Tcl_ObjCmdProc *)polygon_method_hexgrid_create,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::hexagon",(Tcl_ObjCmdProc *)polygon_method_hexagon,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::info",(Tcl_ObjCmdProc *)polygon_method_info,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::intersect",(Tcl_ObjCmdProc *)polygon_method_intersect,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::poly_place",(Tcl_ObjCmdProc *)polygon_method_poly_place,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::rectangle",(Tcl_ObjCmdProc *)polygon_method_rectangle,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::simplify",(Tcl_ObjCmdProc *)polygon_method_simplify,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::segments",(Tcl_ObjCmdProc *)polygon_method_segments,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::squaregrid_create",(Tcl_ObjCmdProc *)polygon_method_squaregrid_create,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::vector_place",(Tcl_ObjCmdProc *)polygon_method_vector_place,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::polygon::within",(Tcl_ObjCmdProc *)polygon_method_within,NULL,NULL);

  Tcl_CreateEnsemble(interp, modPtr->fullName, modPtr, TCL_ENSEMBLE_PREFIX);
  Tcl_Export(interp, modPtr, "[a-z]*", 1);
  
  return TCL_OK;
}

#include "odieInt.h"
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <math.h>

#ifndef IRM_EPSILON
#define IRM_EPSILON 1.19E-07
#endif

#ifndef M_PI
# define M_PI 3.1415926535898
#endif


/*
** LEFT OFF IMPLEMENTING SEGMENT SET AS A TCL OBJECT
** TODO: Write routine to detect a polygon and convert
** to segset and vice versa
**
*/

static Tcl_Interp *local_interp;

#if 0
/*
** Print all vectors in the SegSet.  Used for debugging purposes only.
*/
static void SegPrint(Segment *p, const char *zText){
  printf("%s ", zText);
  if( p ){
    printf("%g,%g -> %g,%g\n", p->from[X_IDX], p->from[Y_IDX], p->to[X_IDX], p->to[Y_IDX]);
  }else{
    printf(" (null)\n");
  }
}
static void SegSetPrint(SegSet *pSet){
  Link *pLink;
  printf("%d vectors:\n", pSet->nSeg);
  for(pLink=pSet->pAll; pLink; pLink=pLink->pNext){
    SegPrint(pLink->pLinkNode, "  ");
  }
}
#endif


/*
** This section implements code for breaking up compartment floorplans
** into triangles.
**
** A floorplan is defined by vectors (x0,y0,x1,y1) which define the
** parameter of each compartment.  The interior of the compartment
** is always to the right of the vector.   Thus the outer boundary
** of the compartment rotates clockwise when viewed from above.
** Compartments may contain holes which are interior voids surrounded 
** by counter-clockwise rotating boundaries.
*/

/*
** Given line segment AB, locate segment BC and return a Vector to
** it.  If there is not BC return NULL.  If there is more than one
** BC return the one that minimizes the angle ABC.
*/
static Segment *SegSetNext(SegSet *pSet, Segment *pAB){
  Link *pX;
  Segment *pBest = 0;
  double angle, bestAngle;
  int cnt = 0;
  int h;
  h = hashVectorXY(pAB->to);

  for(pX=pSet->hashFrom[h]; pX; pX=pX->pNext){
    Segment *pSeg = pX->pLinkNode;
    if( !sameVectorXY(pSeg->from,pAB->to) ) continue;
    /* if(pAB->isBoundary > 1 && pSeg->isBoundary!=0 && pSeg->isBoundary!=pAB->isBoundary) continue; */
    if( cnt==0 ){
      pBest = pSeg;
      bestAngle = fabs(VectorXY_angleOf(pAB->from, pAB->to, pBest->to));
    }else{
      angle = fabs(VectorXY_angleOf(pAB->from, pAB->to, pSeg->to));
      if( angle<bestAngle ){
        bestAngle = angle;
        pBest = pSeg;
      }
    }
    cnt++;
  }
  return pBest;
}

/*
** Remove all segments whose length is less than minLength.  For
** each segment removed, coalesce the inputs into a single new
** VectorXY at the center of the segment.
*/
static void Segset_Bisect_Edges(SegSet *pSet, double minLength){
  Link *pLoop, *pNext;
  double minLen2;

  /* Step 1 - Reduce short edges */
  minLen2 = minLength*minLength;
  
  
  for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
    Segment *p;
    VectorXY from, to, center;
    double vec_dist_sq;

    p = pLoop->pLinkNode;
    pNext = pLoop->pNext;
    if(p->ignore) continue;
    p->ignore=1;
    vec_dist_sq = VectorXY_distance_squared(p->from, p->to);
    VectorXY_Set(from,p->from);
    VectorXY_Set(to,p->to);
    center[X_IDX] = rint(0.5*(from[X_IDX] + to[X_IDX]));
    center[Y_IDX] = rint(0.5*(from[Y_IDX] + to[Y_IDX]));
    
    if(vec_dist_sq>(minLen2*3)) {
      Segment *q;
      /* BISECT */
      p->to[X_IDX]=center[X_IDX];
      p->to[Y_IDX]=center[Y_IDX];
      q=SegSetInsert(pSet,center,to,p->isBoundary);
      if(q) {
        q->midpoint=1;
        q->ignore=1;
      }
    } else if(vec_dist_sq<minLen2 ){
      /* SHRINK */
      SegSetRemove(pSet, p);
      Odie_Free((char *)p);
      for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
        pNext = pLoop->pNext;
        p = pLoop->pLinkNode;
        if( sameVectorXY(p->to,from) || sameVectorXY(p->to,to) ) {
          VectorXY_Set(p->to,center);
        }
        if( sameVectorXY(p->from,from) || sameVectorXY(p->from,to) ) {
          VectorXY_Set(p->from,center);
          SegRelink(pSet, p);
        }
        if( sameVectorXY(p->from,p->to) ) {
          SegSetRemove(pSet, p);
          Odie_Free((char *)p);
        }
      }
      pNext = pSet->pAll;
    }
  }
}

int Segset_Insert_Vectors(Tcl_Interp *interp,SegSet *pSet,int fill,int listLen,Tcl_Obj **listObjPtrs) {
  VECTORXY A,B;
  int i;
  /* Import a flat list, every 4 coordinates x0 y0 x1 y1*/
  if(listLen % 4) {
    Tcl_AppendResult(interp, "Could not interpret coordinates", 0);
    return TCL_ERROR;
  }
  for(i=0;i<listLen;i+=4) {      
    if(Tcl_GetDoubleFromObj(interp, listObjPtrs[i], &A[X_IDX])) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp, listObjPtrs[i+1], &A[Y_IDX])) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp, listObjPtrs[i+2], &B[X_IDX])) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp, listObjPtrs[i+3], &B[Y_IDX])) return TCL_ERROR;
    SegSetInsert(pSet,A,B,fill);
  }
  return TCL_OK;
}

/*
** Internal Utilities
*/

/*
** Routines for handling segments
*/

/*
** Add a new segment to the set
*/
CTHULHU_INLINE Segment *SegSetInsert(
  SegSet *pSet,
  VECTORXY from,
  VECTORXY to,
  u8 isBoundary
){
  int h;
  VectorXY_Round(from);
  VectorXY_Round(to);

  if( sameVectorXY(from,to) ){
    return NULL;
  }
  Segment *p = (Segment *)Odie_Alloc( sizeof(*p) );
  if( p==0 ) {
    return NULL;
  }
  VectorXY_Set(p->from,from);
  VectorXY_Set(p->to,to);

  p->isBoundary = isBoundary;
  p->notOblique = 0;
  LinkInit(p->pAll, p);
  LinkInit(p->pFrom, p);
  LinkInit(p->pSet, p);
  
  LinkInsert(&pSet->pAll, &p->pAll);
  h = hashPoint(p->from);
  LinkInsert(&pSet->hashFrom[h], &p->pFrom);

  pSet->nSeg++;
  pSet->pCurrent = p;
  return p;
}

/*
** Remove a segment from the segment set
*/
CTHULHU_INLINE void SegSetRemove(SegSet *pSet, Segment *p){
  LinkRemove(&p->pAll);
  LinkRemove(&p->pFrom);
  pSet->nSeg--;
  if( pSet->pCurrent==p ){
    pSet->pCurrent = pSet->pAll ? pSet->pAll->pLinkNode : 0;
  }
}

/*
** Call this routine to relink into a segment when the
** Seg.from vector changes.
*/
CTHULHU_INLINE void SegRelink(SegSet *pSet, Segment *p){
  int h;
  LinkRemove(&p->pFrom);
  h = hashPoint(p->from);
  LinkInsert(&pSet->hashFrom[h], &p->pFrom);
}

/*
** Remove all segments from a segment set
*/
CTHULHU_INLINE void SegSetClear(SegSet *pSet){
  while( pSet->pAll ){
    Segment *p;
    assert( pSet->nSeg>0 );
    p=pSet->pAll->pLinkNode;
    SegSetRemove(pSet, p);
    Odie_Free((char *)p);
  }
  assert( pSet->nSeg==0 );
}

/*
** Advance the pSet->pAll pointer so that it is pointing to a different
** segment.
*/
CTHULHU_INLINE void SegSetStep(SegSet *pSet){
  if( pSet->pCurrent ){
    Link *pNext = pSet->pCurrent->pAll.pNext;
    pSet->pCurrent = pNext ? pNext->pLinkNode : 0;
  }
  if( pSet->pCurrent==0 ){
    pSet->pCurrent = pSet->pAll ? pSet->pAll->pLinkNode : 0;
  }
}

static int  segset_method_create (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  local_interp=interp;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N");
    return TCL_ERROR;
  }
  int isnew;
  SegSet *p;
  if( Odie_GetSegmentSetFromObj(interp, objv[1], &p, &isnew) ) return TCL_ERROR;
  if(isnew) {
    Tcl_SetObjResult(interp, Odie_NewSegmentSetObj(p));
  } else {
    Tcl_SetObjResult(interp, objv[1]);
  }
  return TCL_OK;
}


static int  segset_method_add (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  local_interp=interp;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N X0 Y0 X1 Y1");
    return TCL_ERROR;
  }
  int i,j;
  SegSet *pSet;
  VECTORXY A,B;
  double x[4];
  
  if( Odie_GetSegmentGetFromVar(interp, objv[1], &pSet) ) return TCL_ERROR;
  for(i=2;i<objc;i+=4) {
    Segment *found;
    if((objc-i)<4) break;
    for(j=0; j<4; j++){
      if( Tcl_GetDoubleFromObj(interp, objv[i+j], &x[j]) ){
        goto createfail;
      }
    }
    A[X_IDX]=x[0];
    A[Y_IDX]=x[1];
    B[X_IDX]=x[2];
    B[Y_IDX]=x[3];
    /*
    ** Do not insert a vector into the wallset if it
    ** matches a vector already given. It's either redundent
    ** or the edge of a hole
    */
    found=SegSetFind(pSet,A,B);
    if(found) {
      if(found->isBoundary<1) {
        found->isBoundary=1;
      }
    } else {
      SegSetInsert(pSet, A, B, 1);
    }
  }

  Tcl_Obj *objPtr=Odie_NewSegmentSetObj(pSet);
  Tcl_ObjSetVar2(interp,objv[1],NULL,objPtr,0);
  Tcl_SetObjResult(interp, objPtr);

  return TCL_OK;

createfail:
  if(pSet) {
    SegSetClear(pSet);
  }

  return TCL_ERROR;
}

static int  segset_method_subtract (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  local_interp=interp;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "N X0 Y0 X1 Y1");
    return TCL_ERROR;
  }
  int i,j;
  SegSet *pSet=NULL;
  VECTORXY A,B;
  double x[4];

  if( Odie_GetSegmentGetFromVar(interp, objv[1], &pSet) ) return TCL_ERROR;
  for(i=2;i<objc;i+=4) {
    Segment *found;
    if((objc-i)<4) break;
    for(j=0; j<4; j++){
      Tcl_Obj *pObj;
      if( Tcl_GetDoubleFromObj(interp, objv[i+j], &x[j]) ){
        goto createfail;
      }
    }
    A[X_IDX]=x[0];
    A[Y_IDX]=x[1];
    B[X_IDX]=x[2];
    B[Y_IDX]=x[3];
    /*
    ** Do not insert a vector into the wallset if it
    ** matches a vector already given. It's either redundent
    ** or the edge of a hole
    */
    found=SegSetFind(pSet,A,B);
    if(found) {
      if(found->isBoundary<3) {
        found->isBoundary=3;
      }
    } else {
      SegSetInsert(pSet, A, B, 3);
    }
  }
  
  /*
  ** Meander through and clip all segments of isBoundary=1 that touch
  ** an isBoundary=3
  */
  
  Tcl_Obj *objPtr=Odie_NewSegmentSetObj(pSet);
  Tcl_ObjSetVar2(interp,objv[1],NULL,objPtr,0);
  Tcl_SetObjResult(interp, objPtr);
  return TCL_OK;

createfail:
  if(pSet) {
    SegSetClear(pSet);
  }

  return TCL_ERROR;
}


static int  segset_method_difference (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
#ifdef NEVER
  local_interp=interp;
  if( objc!=3 ){
    Tcl_WrongNumArgs(interp, 1, objv, "SEG_POSITIVE SEG_NEGATIVE");
    return TCL_ERROR;
  }
  int created;
  int i,j;
  SegSet *pSet=NULL;
  VECTORXY A,B;
  double x[4];

  Tcl_Obj **varv;
  int varc;
  if(Odie_GetSegmentSetFromObj(interp,objv[1],&pSet,&created)) return TCL_ERROR;

  if (Tcl_ListObjGetElements(interp, objv[2], &varc, &varv) != TCL_OK) {
    goto createfail;
  }
  if(Segset_Insert_Vectors(interp,pSet,3,varc,varv)) {
    goto createfail;
  }
  Link *pLoop, *pNext;
  Link *qLoop, *qNext;

  for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
    Segment *pAB;
    pAB = pLoop->pLinkNode;
    pNext = pLoop->pNext;
    if(pAB->isBoundary>1) continue;

    for(qLoop=pSet->pAll; qLoop; qLoop=pNext){
      Segment *pCD;
      int incident=0;
      pCD = qLoop->pLinkNode;
      qNext = qLoop->pNext;
      if(pCD->isBoundary<3) continue;
      incident=ODIE_Math_LineLineCoincident(
        pAB->from[X_IDX],pAB->from[Y_IDX],
        pAB->to[X_IDX],pAB->to[Y_IDX],
        pCD->from[X_IDX],pCD->from[Y_IDX],
        pCD->to[X_IDX],pCD->to[Y_IDX]
      );
      switch(incident) {
        case 0: continue; /* No overlap */
        case 3:
          /* pAB fits entirely in the range of pBC*/
          SegSetRemove(pSet,pAB);
          continue;
        case 12:
          /* pCD fits entirely in the range of pAB*/
          /* Shorten the first side, and add a stub to represent the other */
          /* THIS CODE IS NOT FINISHED */
          if(VectorXY_distance_squared(pAB->from, pCD->from) < VectorXY_distance_squared(pAB->from, pCD->to)) {
            VectorXY_Set(pAB->from,pCD->from);
            SegRelink(pSet,pAB);
          } else {
            VectorXY_Set(pAB->to,pCD->to);
            SegSetInsert(pSet,pCD->from,pAB->to,1);
          }
          break;
        case 1: {
          /* A is along CD */
          if(VectorXY_distance_squared(pAB->from, pCD->from) < VectorXY_distance_squared(pAB->from, pCD->to)) {
            VectorXY_Set(pAB->from,pCD->from);
          } else {
            VectorXY_Set(pAB->from,pCD->to);
          }
          SegRelink(pSet,pAB);
          break;
        }
        case 2: {
          /* B is along CD */
          /* A is along CD */
          if(VectorXY_distance_squared(pAB->to, pCD->from) < VectorXY_distance_squared(pAB->to, pCD->to)) {
            VectorXY_Set(pAB->to,pCD->from);
          } else {
            VectorXY_Set(pAB->to,pCD->to);
          }
          break;
        }
        case 4: {
          /* C is along AB */
          if(VectorXY_distance_squared(pCD->from, pAB->from) < VectorXY_distance_squared(pCD->from, pAB->to)) {
            VectorXY_Set(pAB->from,pCD->from);
            SegRelink(pSet,pAB);
          } else {
            VectorXY_Set(pAB->to,pCD->from);
          }
          break;
        }
        case 8: {
          /* D is along AB */
          if(VectorXY_distance_squared(pCD->to, pAB->from) < VectorXY_distance_squared(pCD->to, pAB->to)) {
            VectorXY_Set(pAB->from,pCD->to);
            SegRelink(pSet,pAB);
          } else {
            VectorXY_Set(pAB->to,pCD->to);
          }
          break;
        }
      }
    }
  }
  
  
  if(isnew) {
    Tcl_SetObjResult(interp, Odie_NewSegmentSetObj(pSet));
  } else {
    Tcl_InvalidateStringRep(objv[1]);
    Tcl_SetObjResult(interp, objv[1]);
  }  
  return TCL_OK;

createfail:
  if(pSet) {
    SegSetClear(pSet);
  }

  return TCL_ERROR;
#else
  return TCL_OK;
#endif
}

static int  segset_method_rectangle (
  ClientData *dummy,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  double cx, cy, radx,rady;
  Tcl_Obj *pResult=Tcl_NewObj();
  local_interp=interp;

  if( objc != 5 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  radx=radx/2.0;
  rady=rady/2.0;
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}


static int WalkSegments(SegSet *pSet) {
  Link *pLoop, *pNext;
  int changed=0;
  for(pLoop=pSet->pAll; pLoop; pLoop=pNext){
    Segment *pAB, *pBC;
    Link *pRight, *pL2;
    int c;
    /* Find an oblique angle ABC */
    pAB = pLoop->pLinkNode;
    pNext = pLoop->pNext;

    /*
    ** If we are at an oblique for a non boundary
    ** segment, continue
    */
    if( pAB->notOblique ) continue;
    pBC = SegSetNext(pSet, pAB);
    if( pBC==0 ) {
      /*
      ** Remove an orphan wall
      */
      SegSetRemove(pSet, pAB);
      Odie_Free((char *)pAB);
      changed=1;
      continue;
    }

    if( (c = VectorXY_rightOf(pAB->from, pAB->to, pBC->to))>=0 ){
      if( c>0 || !sameVectorXY(pAB->from,pBC->to) ){
        pAB->notOblique = 1;
        continue;
      }
    }

    /* If we reach here, it means that ABC is an oblique angle.
    ** Locate all vertices to the right of AB.
    */
    pRight = 0;
    for(pL2=pSet->pAll; pL2; pL2=pL2->pNext){
      Segment *pX = pL2->pLinkNode;
      if( VectorXY_strictlyRightOf(pAB->from, pAB->to, pX->from)<0 ) continue;
      if( sameVectorXY(pAB->to,pX->from) ) continue;
      pX->score = VectorXY_distance_squared(pAB->to, pX->from);
      pX->isRight = VectorXY_rightOf(pBC->from, pBC->to, pX->from);
      LinkInit(pX->pSet, pX);
      LinkInsert(&pRight, &pX->pSet);
    }
    if( pRight==0 ){
      return TCL_ERROR;
    }
 
    /* pRight is a list of vertices to the right of AB.  Find the
    ** closest vertex X on this list where the line BX does not intersect
    ** any other segment in the polygon.  Then add segments BX and XB.
    */
    while( pRight ){
      Link *pBest=NULL;
      double bestScore;
      int bestRight;
      Segment *pThis,*pX, *pQ;


      /* Search for the "best" vertex.   The best vertex is the
      ** one that is closest.  Though if the vertex is to the left
      ** of BC (and thus would create another oblique angle) then
      ** artificially reduce its score because we would prefer not
      ** to use it.
      */
      pBest = pRight;
      pThis=pBest->pLinkNode;
      bestScore = pThis->score;
      bestRight = pThis->isRight;
      for(pL2=pBest->pNext; pL2; pL2=pL2->pNext){
        int better=0;
        pX = pL2->pLinkNode;
        if( pX->isRight>0 && bestRight <=0 ) {
          better=1;
        } else if ( pX->isRight<=0 && bestRight>0 ) {
          better=0;
        } else if( pX->score<bestScore ){
          better=1;
        }
        if(better) {
          bestScore = pX->score;
          bestRight = pX->isRight;
          pBest = pL2;
        }
      }
      
      

      /* The best vertex is pX */
      pX = pBest->pLinkNode;
      LinkRemove(pBest);

      /* Check to see if BX intersects any segment.  If it does, then
      ** go back and search for a different X
      */
      for(pL2=pSet->pAll; pL2; pL2=pL2->pNext){
        pQ = pL2->pLinkNode;
        if( pQ!=pAB && pQ!=pX
            && VectorXY_intersect(pAB->to, pX->from, pQ->from, pQ->to) ){
          break;
        }
      }
      if( pL2 ) continue;

      /* It did not intersect.  So add BX and XB to the pSet->
      */
      SegSetInsert(pSet, pAB->to, pX->from, 0);
      SegSetInsert(pSet, pX->from, pAB->to, 0);
      pRight = 0;
    }
    changed=1;
    if(!pAB->isBoundary) {
      pNext = pSet->pAll;
    }
  }
  if(changed) {
    return TCL_CONTINUE;
  }
  return TCL_OK;
}


/*
** tclcmd:   convex_subpolygons VECTORS ?MINLENGTH? ?HOLE? ?HOLE? ...
**
** VECTORS is a list of floating-VectorXY values.  Each group of four values
** forms a vector X0,Y0[X_IDX]1,Y1.  The vectors are in no particular order,
** but together they form one or more loops.  Space to the right of each
** vector is within the loop and space to the left is outside.
**
** Loops can be nested.  The outer boundary is formed by a clockwise loop
** of vectors.  Interior holes are formed by counter-clockwise loops.
**
** The output is a list polygons.  Each polygon is a list of 3 or more
** X,Y coordinate pairs.  All polygons are convex and disjoint and they
** together cover the input polygon.
**
** Optionally, the user can specify a series of polygons to be subtracted
** from the main polygon. These are given as an XY list suitable for
** producing a polygon on the tkcanvas
*/
TCL_COMMAND int  segset_method_decompose (
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Tcl_Obj *pSub;   /* A sublist for a single polygon */
  int i, idx, cnt, created;
  SegSet *set;
  double minLen = 0.0;
  local_interp=interp;

  if( objc!=2 && objc!=3 && objc<4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "VECTORS ?MINLENGTH? ?HOLE? ?HOLE? ...");
    return TCL_ERROR;
  }
  if(Odie_GetSegmentSetFromObj(interp,objv[1],&set,&created)) return TCL_ERROR;
  if( objc>2 ){
    if( Tcl_GetDoubleFromObj(interp, objv[2], &minLen) ) return TCL_ERROR;
  }
  /*
  ** Insert the polygons the user specified as
  ** the shape of the holes first
  */
  for(idx=3;idx<objc;idx++) {
    Poly *p;
    int created;
    if(Odie_GetPolygonFromObj(interp,objv[idx],&p,&created)) goto createrror;
    Segset_Insert_Polygon(set,p,idx);
    if(created) Odie_Free((char *)p);
  }

  if( minLen>0.0 ){
    Segset_Bisect_Edges(set,minLen);
  }
  
  cnt=0;
  i=TCL_CONTINUE;
  while(i==TCL_CONTINUE) {
    i=WalkSegments(set);
    cnt++;
    if(cnt>10) {
      break;
    }
  }
  if(i==TCL_CONTINUE) {
    Tcl_AppendResult(interp, "boundary too complex", 0);
    SegSetClear(set);
    return TCL_ERROR;
  }
  if(i==TCL_ERROR) {
    Tcl_AppendResult(interp, "boundary does not enclose a finite space", 0);
    SegSetClear(set);
    return TCL_ERROR;
  }

  /* Now all polygons should be convex.  We just have to generate them. */
  int obtuseangles=0;
  Tcl_Obj *pOut=Tcl_NewObj(); /* The output list */
  //Odie_trace_printf(interp,"NSEG %d\n",set->nSeg);
  while( set->nSeg ){
    VectorXY start;
    
    Segment *pAB, *pBC;
    int valid = 0;
    int cnt = 0;
    
    pAB = set->pAll->pLinkNode;
    start[X_IDX]=pAB->from[X_IDX];
    start[Y_IDX]=pAB->from[Y_IDX];

    /*
    ** Walk along the wallsets, filter out
    ** any that do not include one of the
    ** vectors given as an input of the first
    ** argument
    */
    pSub = Tcl_NewObj();
    while( pAB ){
      pBC = SegSetNext(set, pAB);
      if(pAB->isBoundary < 2) valid=1;
      cnt++;
      if( minLen>=1.0 ){
        Tcl_ListObjAppendElement(0, pSub, Tcl_NewIntObj(pAB->to[X_IDX]));
        Tcl_ListObjAppendElement(0, pSub, Tcl_NewIntObj(pAB->to[Y_IDX]));
      } else {
        Tcl_ListObjAppendElement(0, pSub, Tcl_NewDoubleObj(pAB->to[X_IDX]));
        Tcl_ListObjAppendElement(0, pSub, Tcl_NewDoubleObj(pAB->to[Y_IDX]));
      }
      SegSetRemove(set, pAB);
      if( sameVectorXY(pAB->to,start) ) {
        break;
      }
      Odie_Free((char *)pAB);
      pAB = pBC;
    }
    if( pAB==0 || cnt<3 || !valid){
      Tcl_DecrRefCount(pSub);
    }else{
      Tcl_ListObjAppendElement(0, pOut, pSub);
      //Odie_trace_printf(local_interp,"NEWPOLY %s\n",Tcl_GetString(pSub));
    }
  }

  if(created) {
    SegSetClear(set);
    Odie_Free((char *)set);
  }

  Tcl_SetObjResult(interp, pOut);
  return TCL_OK;

createrror:
  SegSetClear(set);
  return TCL_ERROR;
}

int Odie_Segset_Init(Tcl_Interp *interp) {
  Tcl_Namespace *modPtr;
  
  modPtr=Tcl_FindNamespace(interp,"segset",NULL,TCL_NAMESPACE_ONLY);
  if(!modPtr) {
    modPtr = Tcl_CreateNamespace(interp, "segset", NULL, NULL);
  }
  
  Tcl_CreateObjCommand(interp,"::segset::add",(Tcl_ObjCmdProc *)segset_method_add,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::segset::create",(Tcl_ObjCmdProc *)segset_method_create,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::segset::subtract",(Tcl_ObjCmdProc *)segset_method_subtract,NULL,NULL);
  //Tcl_CreateObjCommand(interp,"::segset::vectors",(Tcl_ObjCmdProc *)segset_method_vectors,NULL,NULL);


  Tcl_CreateObjCommand(interp,"::segset::decompose",(Tcl_ObjCmdProc *)segset_method_decompose,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::segset::rectangle",(Tcl_ObjCmdProc *)segset_method_rectangle,NULL,NULL);

  Tcl_CreateEnsemble(interp, modPtr->fullName, modPtr, TCL_ENSEMBLE_PREFIX);
  Tcl_Export(interp, modPtr, "[a-z]*", 1);
  
  return TCL_OK;
}

/*
** This file is machine generated. Changes will
** be overwritten on the next run of cstruct.tcl
*/
#include "odieInt.h"

/*
** Functions provided by the template
*/


static int  shapes_method_corners (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  double cx, cy, radx,rady;
  
  if( objc != 5 && objc != 9 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy ?x0var y0var x1var y1var?");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  if (objc == 5) {
    Tcl_Obj *pResult=Tcl_NewObj();
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));
    Tcl_SetObjResult(interp, pResult);
    return TCL_OK;
  }
  /*
    Replaces
    set x0 [expr {$cx-$d}]
    set y0 [expr {$cy-$d}]
    set x1 [expr {$cx+$d}]
    set y1 [expr {$cy+$d}]
  */
  
  Tcl_ObjSetVar2(interp,objv[5],NULL,Tcl_NewDoubleObj(cx+radx),0);
  Tcl_ObjSetVar2(interp,objv[6],NULL,Tcl_NewDoubleObj(cy-rady),0);
  Tcl_ObjSetVar2(interp,objv[7],NULL,Tcl_NewDoubleObj(cx-radx),0);
  Tcl_ObjSetVar2(interp,objv[8],NULL,Tcl_NewDoubleObj(cy+rady),0);
  
  return TCL_OK; 
}

static int  shapes_method_drawobj_orientation (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  Tcl_Obj *temp;
  int len;
  double nx=100,ny=0;
  if( objc !=4 ){
    Tcl_WrongNumArgs(interp, 1, objv, "orientation nxvar nyvar");
    return TCL_ERROR;
  }
  
  if(Tcl_ListObjLength(interp,objv[1],&len)) return TCL_ERROR;
  if(len>0) {
    if(Tcl_ListObjIndex(interp, objv[1], 0, &temp)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,temp,&nx)) return TCL_ERROR;
  }
  if(len>1) {
    if(Tcl_ListObjIndex(interp, objv[1], 1, &temp)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,temp,&ny)) return TCL_ERROR;
  }
  Tcl_ObjSetVar2(interp,objv[2],NULL,Tcl_NewDoubleObj(nx),0);
  Tcl_ObjSetVar2(interp,objv[3],NULL,Tcl_NewDoubleObj(ny),0);
  return TCL_OK;
}


static int  shapes_method_poly_hex (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  int i,flip=0;
  double cx, cy, radx,rady;

  Tcl_Obj *pResult=Tcl_NewObj();
  double coords[7][2]= {
    {1.00, 0.00} , {0.50, M_SQRT3_2} , 
    {-0.50, M_SQRT3_2} , {-1.00, -0.00} , 
    {-0.50, -M_SQRT3_2} , {0.50, -M_SQRT3_2},
    {1.00, 0.00} 
  };
  if( objc != 5 && objc != 6){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy ?flip?");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  if(objc==6) {
    if(Tcl_GetBooleanFromObj(interp,objv[5],&flip)) return TCL_ERROR;
  }
  radx=radx/2.0;
  rady=rady/2.0;

  for(i=0;i<6;i++) {
    if(flip) {
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i+1][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i+1][0]));
    } else {
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i][1]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx*coords[i+1][0]));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady*coords[i+1][1]));
    }
  }

  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  shapes_method_poly_place (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  /*
  ** Apply Matrices
  */
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double zoom;
  double matA[6] = {1.0,0.0,0.0,1.0,0.0,0.0};
  double centerx,centery,normalx,normaly,angle;

  if( objc < 8 ){
      Tcl_WrongNumArgs(interp, 1, objv, "zoom centerx centery normalx normaly x1 y1 ?x2 y2?...");
      return TCL_ERROR;
  }   
  if(Tcl_GetDoubleFromObj(interp,objv[1],&zoom)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&centerx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&centery)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&normalx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[5],&normaly)) return TCL_ERROR;

  angle=atan2(normaly,normalx);
  matA[0]=cos(angle);
  matA[1]=sin(angle);
  matA[2]=-sin(angle);
  matA[3]=cos(angle);
  matA[4]=0.0;
  matA[5]=0.0;
  double startx,starty,prevx,prevy;

  i=6;
  {
    double x,y,sx,sy,newx,newy;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
    
    sx=(x/zoom);
    sy=(y/zoom);
    newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
    newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

    startx=newx;
    starty=newy;
    prevx=newx;
    prevy=newy;
  }
  
  for(i=8;i<objc;i+=2) {
    double x,y,sx,sy,newx,newy;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
    
    sx=(x/zoom);
    sy=(y/zoom);
    newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
    newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevy));  
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newy));       

    prevx=newx;
    prevy=newy;
  }
  if(startx != prevx && starty!= prevy) {
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(prevy));  
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(startx));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(starty));  
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  shapes_method_polygon_to_vectors (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double px,py,fx,fy;

  if(objc<6 || (objc-1)%2!=0) {
      Tcl_AppendResult(interp, "arguments should contain at least 6 and "
           " a multiple of 2 values", 0);
      return TCL_ERROR;

  }
  if(Tcl_GetDoubleFromObj(interp,objv[1],&px)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&py)) return TCL_ERROR;
  fx=px;
  fy=py;
  for(i=3;i<objc;i+=2) {
    double x,y;
    if(i+1>objc) break;
    if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
    if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(px));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(py));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(x));
    Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(y));
    px=x;
    py=y;
  }
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(px));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(py));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(fx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(fy));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  shapes_method_rectangle_as_polygon (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  double cx, cy, radx,rady;
  Tcl_Obj *pResult=Tcl_NewObj();

  if( objc != 5 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  radx=radx/2.0;
  rady=rady/2.0;
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  shapes_method_rectangle_as_vectors (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {


  double cx, cy, radx,rady;
  Tcl_Obj *pResult=Tcl_NewObj();

  if( objc != 5 ){
    Tcl_WrongNumArgs(interp, 1, objv, "cx cy dimx dimy");
    return TCL_ERROR;
  }
  
  if(Tcl_GetDoubleFromObj(interp,objv[1],&cx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&cy)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&radx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&rady)) return TCL_ERROR;
  radx=radx/2.0;
  rady=rady/2.0;
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx+radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));
  
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy+rady));

  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cx-radx));
  Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(cy-rady));
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

static int  shapes_method_vector_place (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
  /*
  ** Apply Matrices
  */
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double zoom;
  double matA[6] = {1.0,0.0,0.0,1.0,0.0,0.0};
  double centerx,centery,normalx,normaly,angle;

  if( objc < 8 ){
      Tcl_WrongNumArgs(interp, 1, objv, "zoom centerx centery normalx normaly x1 y1 ?x2 y2?...");
      return TCL_ERROR;
  }   
  if(Tcl_GetDoubleFromObj(interp,objv[1],&zoom)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&centerx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&centery)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[4],&normalx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[5],&normaly)) return TCL_ERROR;

  angle=atan2(normaly,normalx);
  matA[0]=cos(angle);
  matA[1]=sin(angle);
  matA[2]=-sin(angle);
  matA[3]=cos(angle);
  matA[4]=0.0;
  matA[5]=0.0;

  
  for(i=6;i<objc;i+=2) {
      double x,y,sx,sy,newx,newy;
      if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
      if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
      
      sx=(x/zoom);
      sy=(y/zoom);
      newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
      newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newx));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newy));
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}


static int  shapes_method_canvas (
  ClientData *simulator,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
) {
   /*
  ** Apply Matrices
  */
  Tcl_Obj *pResult=Tcl_NewObj();
  int i;
  double zoom;
  double matA[6] = {1.0,0.0,0.0,1.0,0.0,0.0};
  double centerx,centery,width,height,angle=0.0;

  if( objc < 6 || objc>7 ){
      Tcl_WrongNumArgs(interp, 1, objv, "zoom centerx centery ?angle?");
      return TCL_ERROR;
  }   
  if(Tcl_GetDoubleFromObj(interp,objv[1],&zoom)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&centerx)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&centery)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[2],&width)) return TCL_ERROR;
  if(Tcl_GetDoubleFromObj(interp,objv[3],&height)) return TCL_ERROR;
  if(objc==7) {
    if(Tcl_GetDoubleFromObj(interp,objv[4],&angle)) return TCL_ERROR;
    
  }
  matA[0]=cos(angle);
  matA[1]=sin(angle);
  matA[2]=-sin(angle);
  matA[3]=cos(angle);
  matA[4]=0.0;
  matA[5]=0.0;

  
  for(i=6;i<objc;i+=2) {
      double x,y,sx,sy,newx,newy;
      if(Tcl_GetDoubleFromObj(interp,objv[i],&x)) return TCL_ERROR;
      if(Tcl_GetDoubleFromObj(interp,objv[i+1],&y)) return TCL_ERROR;
      
      sx=(x/zoom);
      sy=(y/zoom);
      newx=matA[0]*sx+matA[1]*sy+matA[4]+centerx;
      newy=matA[2]*sx+matA[3]*sy+matA[5]+centery;

      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newx));
      Tcl_ListObjAppendElement(interp,pResult,Tcl_NewDoubleObj(newy));
  }
  Tcl_SetObjResult(interp, pResult);
  return TCL_OK;
}

int Odie_Shapes_Init(Tcl_Interp *interp) {
  Tcl_Namespace *modPtr;

  modPtr=Tcl_FindNamespace(interp,"shapes",NULL,TCL_NAMESPACE_ONLY);
  if(!modPtr) {
    modPtr = Tcl_CreateNamespace(interp, "shapes", NULL, NULL);
  }

  Tcl_CreateObjCommand(interp,"::shapes::corners",(Tcl_ObjCmdProc *)shapes_method_corners,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::drawobj_orientation",(Tcl_ObjCmdProc *)shapes_method_drawobj_orientation,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::poly_hex",(Tcl_ObjCmdProc *)shapes_method_poly_hex,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::poly_place",(Tcl_ObjCmdProc *)shapes_method_poly_place,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::polygon_to_vectors",(Tcl_ObjCmdProc *)shapes_method_polygon_to_vectors,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::rectangle_as_polygon",(Tcl_ObjCmdProc *)shapes_method_rectangle_as_polygon,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::rectangle_as_vectors",(Tcl_ObjCmdProc *)shapes_method_rectangle_as_vectors,NULL,NULL);
  Tcl_CreateObjCommand(interp,"::shapes::vector_place",(Tcl_ObjCmdProc *)shapes_method_vector_place,NULL,NULL);

  Tcl_CreateEnsemble(interp, modPtr->fullName, modPtr, TCL_ENSEMBLE_PREFIX);
  Tcl_Export(interp, modPtr, "[a-z]*", 1);
  
  return TCL_OK;
}

/*
** This widget translates 3-D coordinates onto a flat canvas by splitting
** the 3-D space into layers and stacking the layers on the canvas.
**
** The layers are decks of the ship.  The highest layer (or deck) is drawn
** at the top of the page.  The next layer down is drawn below the top layer.
** and so forth down the canvas.  In other words, the 3D object is drawn
** by showing a set of 2D slices where each slice is viewed from above.
**
** The original 3D coordinates are called "actual" coordinates.  When
** translated into the 2D canvas they are called "canvas" coordinates.
**
** The actual coordinate system is right-handed.  The X axis increases to
** the right.  The Y axis increases going up.  The Z axis comes out of the
** page at the viewer.  The canvas coordinate system is left-handed.  The
** X axis increase to the right but the Y axis increases going down.
**
** A slicer is a object with methods.  The details of the available
** methods and what each does are described in comments before the
** implementation of each method.
*/
#include "odieInt.h"
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h>

/*
** This routine is called when a slicer is deleted.  All the memory and
** other resources allocated by this slicer is recovered.
*/
static void destroySlicer(void *pArg){
  Slicer *p = (Slicer*)pArg;
  int i;
  for(i=0; i<p->nSlice; i++){
    Odie_Free((char *)p->a[i].zName);
    Odie_Free((char *)p->a[i].xz);
  }
  Odie_Free((char *)p->a);
  Odie_Free((char *)p);
}

static int Location_FromTclObj(Tcl_Interp *interp, Tcl_Obj *pList,int *did,double *x,double *y) {
  int listlen;
  Tcl_Obj **elist;
  double z;
  if(Tcl_ListObjGetElements(interp,pList,&listlen,&elist)) {
    return TCL_ERROR;
  }
  if(listlen < 3 || listlen > 4) {
    Tcl_AppendResult(interp, "Could not interpret location ", Tcl_GetString(pList), 0);
    return TCL_ERROR;
  }
  if( Tcl_GetIntFromObj(interp, elist[0], did) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, elist[1], x) ) return TCL_ERROR;
  if( Tcl_GetDoubleFromObj(interp, elist[2], y) ) return TCL_ERROR;
  return TCL_OK;
}


static double xCanvasToActual(Slicer *p,struct OneSlice *pS, double cx){
   double ax;
  if(p->nSlice==1) {
    ax = cx*p->rZoom;
  } else {
    ax = cx*p->rZoom - pS->rXShift;
  }   
  return ax;   
}

static double yCanvasToActual(Slicer *p,struct OneSlice *pS, double cy){
  double ay;
  if(p->nSlice==1) {
    ay=-cy*p->rZoom;
  } else {
    ay = pS->mxY + (pS->mnY-pS->mxY)*(cy-pS->top)/(pS->btm-pS->top);
  }
  return ay;
}

/*
** Convert a Y coordinate from actual to canvas coordinates for a
** given deck.
*/
static double xActualToCanvas(Slicer *p,struct OneSlice *pS, double ax){
  double cx;
  if(p->nSlice==1) {
    cx = ax/p->rZoom;
  } else {
    cx = (ax+pS->rXShift)/p->rZoom;
  }
  return cx;
}

/*
** Convert a Y coordinate from actual to canvas coordinates for a
** given deck.
*/
static double yActualToCanvas(Slicer *p,struct OneSlice *pS, double ay){
  double cy;
  if(p->nSlice==1) {
    cy=-ay/p->rZoom;
  } else {
    cy = pS->top + (pS->btm-pS->top)*(ay-pS->mxY)/(pS->mnY-pS->mxY);
  }
  return cy;
}

/*
** Return the height above baseline for the deck location rX.
*/
static double deckHeight(struct OneSlice *p, double rX){
  int i;
  double *xz;
  if(!p) {
    return 0.0;
  }
  if( p->nXZ<4 ){
    return p->z;
  }
  xz = p->xz;
  if( rX<=xz[0] ){
    return xz[1];
  }
  for(i=2; i<p->nXZ; i+=2){
    if( rX<=xz[i] ){
      assert( xz[i]>xz[i-2] );
      return xz[i-1] + (xz[i+1]-xz[i-1])*(rX-xz[i-2])/(xz[i]-xz[i-2]);
    }
  }
  return xz[p->nXZ-1];
}

/*
** Return a pointer to the particular deck at actual coordinates
** X, Z
*/
static struct OneSlice *deckAt(Slicer *p, double rX, double rZ){
  int i;
  struct OneSlice *pBest;
  double bestHeight;
  if( p->nSlice==0 ) return 0;
  pBest = &p->a[0];
  bestHeight = deckHeight(pBest, rX);
  for(i=1; i<p->nSlice; i++){
    double dh = deckHeight(&p->a[i],rX);
    if( dh>bestHeight && dh<=rZ ){
      pBest = &p->a[i];
      bestHeight = dh;
    }
  }
  return pBest;
}


/*
** Return a pointer to the deck immediately above the given deck.
** Return NULL if the given deck is topmost.
*/
static struct OneSlice *deckAbove(Slicer *p, struct OneSlice *pRef){
  int i;
  struct OneSlice *pBest = 0;
  if( p->nSlice==0 ) return 0;
  for(i=0; i<p->nSlice; i++){
    struct OneSlice *pTest = &p->a[i];
    if( pTest->z<=pRef->z ) continue;
    if( pBest==0 || pBest->z>pTest->z ){
      pBest = pTest;
    }
  }
  return pBest;
}
/*
** Return a pointer to the deck immediately above the given deck.
** Return NULL if the given deck is topmost.
*/
static struct OneSlice *deckBelow(Slicer *p, struct OneSlice *pRef){
  int i;
  struct OneSlice *pBest = 0;
  if( p->nSlice==0 ) return 0;
  for(i=0; i<p->nSlice; i++){
    struct OneSlice *pTest = &p->a[i];
    if( pTest->z>=pRef->z ) continue;
    if( pBest==0 || pBest->z<pTest->z ){
      pBest = pTest;
    }
  }
  return pBest;
}


/*
** Recompute the values of p->a[].top and p->a[].btm for all slices in
** the given slicer.
*/
static void computeTopAndBottom(Slicer *p){
  int i;
  double rY = 0.0;
  double rBound = -9.9e99;
  for(i=p->nSlice-1; i>=0; i--){
    double h = (p->a[i].mxY - p->a[i].mnY)/p->rZoom;
    p->a[i].upperbound = rBound;
    p->a[i].top = rY;
    p->a[i].btm = rY + h;
    rY = p->a[i].btm + 0.3*h;
    rBound = p->a[i].btm + 0.15*h;
  }
    
  if( p->nSlice==0 ) return;

  /* Calculate the above and below for each deck */

  for(i=0; i<p->nSlice; i++){
    struct OneSlice *pThis = &p->a[i];
    struct OneSlice *pBest = 0;

    pBest=deckAbove(p,pThis);
    if(pBest) {
      pThis->above=pBest->did;
    } else {
      pThis->above=0;
    }
    pBest=deckBelow(p,pThis);
    if(pBest) {
      pThis->below=pBest->did;
    } else {
      pThis->below=0;
    }
  }
}

/*
** pObj is either the name of a deck or a Z coordinate.  If it is a
** deck name, find the deck and write a pointer to it in *ppS.  If
** it is a Z coordinate, use that coordinate together with rX to 
** find the deck and write it into *ppS.  If an error occurs, put
** an error message on the TCL interpreter and return TCL_ERROR.
** Return TCL_OK on success.
*/
static int getDeck(
  Tcl_Interp *interp,     /* Put error messages here */
  Slicer *p,              /* The slicer */
  double rX,              /* X coord used to find deck if pObj is a Z coord */
  Tcl_Obj *pObj,          /* Either a deck name or a Z coordinate */
  struct OneSlice **ppS   /* Write the slice pointer here */
){
  double rZ;
  const char *zName;
  int i;
  if(p->nSlice==1) {
    *ppS=&p->a[0];
    return TCL_OK;
  }
  if( Tcl_GetDoubleFromObj(0, pObj, &rZ)==TCL_OK ){
    *ppS = deckAt(p, rX, rZ);
    return TCL_OK;
  }
  zName = Tcl_GetStringFromObj(pObj, 0);
  for(i=0; i<p->nSlice; i++){
    if( strcmp(zName, p->a[i].zName)==0 ){
      *ppS = &p->a[i];
      return TCL_OK;
    }
  }
  Tcl_AppendResult(interp, "no such deck: ", zName, 0);
  return TCL_ERROR;
}


/*
** pObj is either the name of a deck or a Z coordinate.  If it is a
** deck name, find the deck and write a pointer to it in *ppS.  If
** it is a Z coordinate, use that coordinate together with rX to 
** find the deck and write it into *ppS.  If an error occurs, put
** an error message on the TCL interpreter and return TCL_ERROR.
** Return TCL_OK on success.
*/
static int getDeckId(
  Tcl_Interp *interp,     /* Put error messages here */
  Slicer *p,              /* The slicer */
  Tcl_Obj *pObj,          /* Either a deck name or a Z coordinate */
  struct OneSlice **ppS   /* Write the slice pointer here */
){
  int did;
  const char *zName;
  int i;
  if(p->nSlice==1) {
    *ppS=&p->a[0];
    return TCL_OK;
  }
  if( Tcl_GetIntFromObj(interp, pObj, &did)==TCL_OK ){
    for(i=0; i<p->nSlice; i++){
      if( did == p->a[i].did ){
        *ppS = &p->a[i];
        return TCL_OK;
      }
    }
    Tcl_AppendResult(interp, "no such deckid: ", 0);
    return TCL_ERROR;
  }
  zName = Tcl_GetStringFromObj(pObj, 0);
  for(i=0; i<p->nSlice; i++){
    if( strcmp(zName, p->a[i].zName)==0 ){
      *ppS = &p->a[i];
      return TCL_OK;
    }
  }
  Tcl_AppendResult(interp, "no such deck: ", zName, 0);
  return TCL_ERROR;
}

static inline struct OneSlice *getDeck_FromInt(
  Slicer *p,
  int did
) {
  int i;
  if(p->nSlice==1) {
    return &p->a[0];
  }
  for(i=0; i<p->nSlice; i++){
    if( did == p->a[i].did ){
      return &p->a[i];
    }
  }
  return NULL;
}

/*
** Methods
*/

static int slicer_drawline_do(
  Tcl_Interp *interp,
  Slicer *p,
  Tcl_Obj *canvas,
  Tcl_Obj *tagname,
  Tcl_Obj *tagname_transdeck,
  Tcl_Obj *tagname_penetration,
  int coord_count,
  int *deckCoord,double *xCoord,double *yCoord,
  struct OneSlice **apDeck
) {
  int i;

  Tcl_Obj *pVTag;             /* The "sNNN" tag added to all line segments */
  const char *zXTag;          /* Trans-deck tag (dashed lines) */
  const char *zPTag;          /* Deck-penetraction tag */
  Tcl_Obj *aLineArg[20];      /* Element of "create line" TCL command */
  int nLineArg;               /* Number of used entries in aLineArg[] */
  Tcl_Obj *aPenArg[20];       /* Cmd to draw deck penetractions */
  int nPenArg;
  zXTag = Tcl_GetStringFromObj(tagname_transdeck, 0);
  zPTag = Tcl_GetStringFromObj(tagname_penetration, 0);
  
  aLineArg[0] = canvas;
  aLineArg[1] = ODIE_CONSTANT_STRING("create");
  aLineArg[2] = ODIE_CONSTANT_STRING("line");
  for(i=3; i<=6; i++){
    aLineArg[i] = Tcl_NewObj();
  }
  aLineArg[7] = ODIE_CONSTANT_STRING("-tags");
  for(i=0; i<=7; i++){
    Tcl_IncrRefCount(aLineArg[i]);
  }
  nLineArg = 9;

  if( zPTag[0]==0 ){
    nPenArg = 0;
  }else{
    aPenArg[0] = canvas;
    aPenArg[1] = ODIE_CONSTANT_STRING("create");
    aPenArg[2] = ODIE_CONSTANT_STRING("oval");
    for(i=3; i<=6; i++){
      aPenArg[i] = Tcl_NewObj();
    }
    aPenArg[7] = ODIE_CONSTANT_STRING("-tags");
    for(i=0; i<=7; i++){
      Tcl_IncrRefCount(aPenArg[i]);
    }
    nPenArg = 9;
  }

  for(i=1; i<coord_count; i++){
    char zBuf[30];
    double x0, y0, x1, y1, Y0, Y1;
    sprintf(zBuf, "s%d", i);
    pVTag = Tcl_NewStringObj(zBuf, -1);
    Tcl_IncrRefCount(pVTag);
    if( apDeck[i]!=apDeck[i-1] ){
      x0 = xActualToCanvas(p,apDeck[i-1],xCoord[i-1]);
      x1 = xActualToCanvas(p,apDeck[i],xCoord[i]);
      
      y0 = yActualToCanvas(p,apDeck[i-1], Y0=yCoord[i-1]);
      y1 = yActualToCanvas(p,apDeck[i], Y1=yCoord[i]);
      if( zXTag[0]!=0 ){
        Tcl_SetDoubleObj(aLineArg[3], x0);
        Tcl_SetDoubleObj(aLineArg[4], y0);
        Tcl_SetDoubleObj(aLineArg[5], x1);
        Tcl_SetDoubleObj(aLineArg[6], y1);
        aLineArg[8] = Tcl_DuplicateObj(tagname_transdeck);
        Tcl_ListObjAppendElement(0, aLineArg[8], pVTag);
        Tcl_IncrRefCount(aLineArg[8]);
        Tcl_EvalObjv(interp, nLineArg, aLineArg, 0);
        Tcl_DecrRefCount(aLineArg[8]);
      }
      continue;
    }
    
    x0 = xActualToCanvas(p,apDeck[i],xCoord[i-1]);
    x1 = xActualToCanvas(p,apDeck[i],xCoord[i]);
    
    y0 = yActualToCanvas(p,apDeck[i], yCoord[i-1]);
    y1 = yActualToCanvas(p,apDeck[i], yCoord[i]);
    Tcl_SetDoubleObj(aLineArg[3], x0);
    Tcl_SetDoubleObj(aLineArg[4], y0);
    Tcl_SetDoubleObj(aLineArg[5], x1);
    Tcl_SetDoubleObj(aLineArg[6], y1);
    aLineArg[8] = Tcl_DuplicateObj(tagname);
    Tcl_ListObjAppendElement(0, aLineArg[8], pVTag);
    Tcl_DecrRefCount(pVTag);
    Tcl_IncrRefCount(aLineArg[8]);
    Tcl_EvalObjv(interp, nLineArg, aLineArg, 0);
    Tcl_DecrRefCount(aLineArg[8]);
  }
  for(i=0; i<=7; i++){
    Tcl_DecrRefCount(aLineArg[i]);
    if( i<nPenArg ) Tcl_DecrRefCount(aPenArg[i]);
  }
  for(i=9; i<nLineArg; i++){
    Tcl_DecrRefCount(aLineArg[i]);
  }
  for(i=9; i<nPenArg; i++){
    Tcl_DecrRefCount(aPenArg[i]);
  }
  return TCL_OK;
}

/*
** tclmethod:  SLICER drawline_dxyz CANVAS START PATH END TAG TRANSDECK-TAG P-TAG
** title:   Return TK to draw a line on a canvas
*/ 
static int slicer_method_drawline(
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Slicer *p = (Slicer*)pArg;

  int i, j=0, totalN, n, mode;
  double *xCoord,*yCoord;             /* Actual coordinates of line to draw */
  int *deckCoord;
  struct OneSlice **apDeck;   /* Array of all decks */
  
  if( objc!=6 ){
    Tcl_WrongNumArgs(interp, 1, objv, 
        "CANVAS PATH TAG TRANSDECK-TAG PENETRATION-TAG");
    return TCL_ERROR;
  }
  if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;

  totalN=n+2;
  xCoord = (double *)Odie_Alloc(sizeof(xCoord[0])*totalN);
  yCoord = (double *)Odie_Alloc(sizeof(yCoord[0])*totalN);
  deckCoord = (int *)Odie_Alloc(sizeof(deckCoord[0])*(totalN));
  apDeck = (struct OneSlice **)Odie_Alloc(sizeof(apDeck[0])*(totalN));
  
  if( xCoord==0 || yCoord == 0 || deckCoord == 0 || apDeck == 0 ) {
    return TCL_ERROR;
  }
  
  j=0;
  for(i=0; i<n; i++){
    Tcl_Obj *deckObj,*xObj,*yObj;
    Tcl_ListObjIndex(0, objv[2], i, &deckObj);
    if( Tcl_GetIntFromObj(interp, deckObj, &deckCoord[j]) || i>(n-4) ) {
      if(Location_FromTclObj(interp,deckObj,&deckCoord[j],&xCoord[j],&yCoord[j])) {
        goto badRoute;
      }
    } else {
      if(i++ >= n) goto badRoute; 
      Tcl_ListObjIndex(0, objv[2], i, &xObj);
      if(i++ >= n) goto badRoute; 
      Tcl_ListObjIndex(0, objv[2], i, &yObj);
      i++;
      if( Tcl_GetDoubleFromObj(interp, xObj, &xCoord[j]) ) goto badRoute;
      if( Tcl_GetDoubleFromObj(interp, yObj, &yCoord[j]) )  goto badRoute;
    }
    if(deckCoord[j]==deckCoord[j-1]) {
      apDeck[j]=apDeck[j-1];
    } else {
      apDeck[j]=getDeck_FromInt(p,deckCoord[j]); if(!apDeck[j]) goto badRoute;
    }
    j++;
  }
  slicer_drawline_do(interp,p,objv[1],objv[3],objv[4],objv[5],j,deckCoord,xCoord,yCoord,apDeck);
  Odie_Free((char *)xCoord);
  Odie_Free((char *)yCoord);
  Odie_Free((char *)deckCoord);
  Odie_Free((char *)apDeck);
  return TCL_OK;

badRoute:
  Odie_Free((char *)xCoord);
  Odie_Free((char *)yCoord);
  Odie_Free((char *)deckCoord);
  Odie_Free((char *)apDeck);
  return TCL_ERROR;
}


/*
** tclmethod:  SLICER drawline_dxyz START PATH END
** title:   Return TK to draw a line on a canvas
*/ 
static int slicer_method_link_coords(
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Slicer *p = (Slicer*)pArg;

  int i, j=0, totalN, n, mode;
  double *xCoord,*yCoord;             /* Actual coordinates of line to draw */
  int *deckCoord;
  struct OneSlice **apDeck;   /* Array of all decks */

  if( objc!=8 ){
    Tcl_WrongNumArgs(interp, 1, objv, 
        "START PATH END");
    return TCL_ERROR;
  }
  if( Tcl_ListObjLength(interp, objv[3], &n) ) return TCL_ERROR;

  totalN=n+2;
  xCoord = (double *)Odie_Alloc(sizeof(xCoord[0])*totalN);
  yCoord = (double *)Odie_Alloc(sizeof(yCoord[0])*totalN);
  deckCoord = (int *)Odie_Alloc(sizeof(deckCoord[0])*(totalN));
  apDeck = (struct OneSlice **)Odie_Alloc(sizeof(apDeck[0])*(totalN));
  
  if( xCoord==0 || yCoord == 0 || deckCoord == 0 || apDeck == 0 ) {
    return TCL_ERROR;
  }
  
  j=0;
  if(Location_FromTclObj(interp,objv[2],&deckCoord[j],&xCoord[j],&yCoord[j])) {
    goto badRoute;
  }
  apDeck[j]=getDeck_FromInt(p,deckCoord[j]); if(!apDeck[j]) goto badRoute;

  j++;

  for(i=0; i<n; i++){
    Tcl_Obj *deckObj,*xObj,*yObj;
    Tcl_ListObjIndex(0, objv[3], i, &deckObj);
    if( Tcl_GetIntFromObj(interp, deckObj, &deckCoord[j]) || i>(n-4) ) {
      if(Location_FromTclObj(interp,deckObj,&deckCoord[j],&xCoord[j],&yCoord[j])) {
        goto badRoute;
      }
    } else {
      Tcl_ListObjIndex(0, objv[3], i+1, &xObj);
      Tcl_ListObjIndex(0, objv[3], i+2, &yObj);
      if( Tcl_GetDoubleFromObj(interp, xObj, &xCoord[j]) ) goto badRoute;
      if( Tcl_GetDoubleFromObj(interp, yObj, &yCoord[j]) )  goto badRoute;
      i+=3;
    }
    if(deckCoord[j]==deckCoord[j-1]) {
      apDeck[j]=apDeck[j-1];
    } else {
      apDeck[j]=getDeck_FromInt(p,deckCoord[j]); if(!apDeck[j]) goto badRoute;
    }
    j++;
  }

  if(Location_FromTclObj(interp,objv[4],&deckCoord[j],&xCoord[j],&yCoord[j])) {
    goto badRoute;
  }
  apDeck[j]=getDeck_FromInt(p,deckCoord[j]); if(!apDeck[j]) goto badRoute;
  j++;
  slicer_drawline_do(interp,p,objv[1],objv[5],objv[6],objv[7],j,deckCoord,xCoord,yCoord,apDeck);
  Odie_Free((char *)xCoord);
  Odie_Free((char *)yCoord);
  Odie_Free((char *)deckCoord);
  Odie_Free((char *)apDeck);
  return TCL_OK;

badRoute:
  Odie_Free((char *)xCoord);
  Odie_Free((char *)yCoord);
  Odie_Free((char *)deckCoord);
  Odie_Free((char *)apDeck);
  return TCL_ERROR;
}

static inline double Location_zabs(Slicer *p,struct OneSlice *pS,double x0,double dheight) {
  if(dheight >= 0) {
    return dheight+deckHeight(pS,x0);
  }
  struct OneSlice *pAbove;
  pAbove=deckAbove(p,pS);
  if(pAbove) {
    return deckHeight(pAbove,x0)+dheight;
  }
  return deckHeight(pS,x0)+p->upper_height+dheight;
}

static inline double Location_zdeck(Slicer *p,struct OneSlice *pS,double x0,double dheight) {
  if(dheight >= 0) {
    return dheight;
  }
  struct OneSlice *pAbove;
  pAbove=deckAbove(p,pS);
  if(pAbove) {
    return deckHeight(pAbove,x0)-deckHeight(pS,x0)+dheight;
  }
  return p->upper_height+dheight;
}

/*
** This routine runs when a method is executed against a slicer
*/
static int slicerMethodProc(
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Slicer *p = (Slicer*)pArg;

#if 0
  /* For debugging....
  ** Print each wallset command before it is executed.
  */
  { int i;
    for(i=0; i<objc; i++){
      printf("%s%c", Tcl_GetStringFromObj(objv[i], 0), i<objc-1 ? ' ' : '\n');
    }
  }
#endif

  /* The following bit of magic implements a switch() statement that
  ** invokes the correct code based on the value of objv[1].  The
  ** mktclopts.tcl script scans this source file looking for "case"
  ** statements then generates the "wallset.h" include file that contains
  ** all of the code necessary to implement the switch.
  **
  ** In this way, we can add new methods to the command simply by adding
  ** new cases within the switch.  All the related switch code is
  ** regenerated automatically.
  */
  {
#include "slicer_cases.h"
  {

  /* tclmethod:  SLICER above NAME
  ** title:   Return the deck above NAME
  */
  case SLICER_ABOVE: {
    struct OneSlice *pDeck;
    if( objc!=3 &&  objc!=5 && objc!=6){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME ?x y? ?zoff?");
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    if(objc==3) {
      Tcl_SetObjResult(interp, Tcl_NewIntObj(pDeck->above));
    } else {
      Tcl_Obj *pResult=Tcl_NewObj();

      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewIntObj(pDeck->above));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[3]));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[4]));
      if(objc==6) {
        Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[5]));
      }
      Tcl_SetObjResult(interp, pResult);
    }
    return TCL_OK;
  }
  
  /* tclmethod:  SLICER below NAME
  ** title:   Return the deck below NAME
  */
  case SLICER_BELOW: {
    struct OneSlice *pDeck;
    if( objc!=3 &&  objc!=5 && objc!=6){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME ?x y? ?zoff?");
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    if(objc==3) {
      Tcl_SetObjResult(interp, Tcl_NewIntObj(pDeck->below));
    } else {
      Tcl_Obj *pResult=Tcl_NewObj();

      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewIntObj(pDeck->below));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[3]));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[4]));
      if(objc==6) {
        Tcl_ListObjAppendElement(interp, pResult, Tcl_DuplicateObj(objv[5]));
      }
      Tcl_SetObjResult(interp, pResult);
    }
    return TCL_OK;
  }
  case SLICER_DECKID_TO_NAME: {
    struct OneSlice *pDeck;
    if( objc!=3 ) {
      Tcl_WrongNumArgs(interp, 2, objv, "ID");
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    Tcl_SetObjResult(interp, ODIE_CONSTANT_STRING(pDeck->zName));
    return TCL_OK;
  }
  case SLICER_DECKNAME_TO_ID: {
    struct OneSlice *pDeck;
    if( objc!=3 ) {
      Tcl_WrongNumArgs(interp, 2, objv, "NAME");
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    Tcl_SetObjResult(interp, Tcl_NewIntObj(pDeck->did));
    return TCL_OK;
  }

  /*
  ** tclmethod:  SLICER xyz_to_location X Y Z
  ** title:   Convert from X Y Z to IRM Coordinates
  **
  ** The ABOVE-DECK parameter, if present, is the height above the deck.
  */  
  case SLICER_XYZ_TO_LOCATION: {
    double x0, y0, z0, dheight;
    Tcl_Obj *pResult;
    struct OneSlice *pS;
    int i;
    if(objc==2) {
      Tcl_SetObjResult(interp, Tcl_NewObj());
      return TCL_OK;
    }
    if( objc < 5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y Z ?x y z?");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    pResult = Tcl_NewObj();

    for(i=2;i<objc;i+=3) {
      if( Tcl_GetDoubleFromObj(interp, objv[i], &x0) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[i+1], &y0) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[i+2], &z0) ) return TCL_ERROR;
      pS = deckAt(p, x0, z0);
      dheight=z0-deckHeight(pS,x0);
    
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pS->did));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(x0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(y0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(dheight));
    }

    Tcl_SetObjResult(interp, pResult);
    break;
  }
  case SLICER_LOCATION_TO_XYZ: {
    int i;
    double x0, y0, z0, dheight;
    Tcl_Obj *pResult;
    struct OneSlice *pS,*pAbove;
    if(objc==2) {
      Tcl_SetObjResult(interp, Tcl_NewObj());
      return TCL_OK;
    }
    if( objc!= 5 && objc < 6 ){
      Tcl_WrongNumArgs(interp, 2, objv, "DECK X Y ?ZOFF? ...");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if(objc<7) {
      if( getDeckId(interp,p,objv[2],&pS) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[4], &y0) ) return TCL_ERROR;
      if(objc==6) {
        if( Tcl_GetDoubleFromObj(NULL, objv[5], &dheight) ) {
          dheight=0.0;
        }
      } else {
        dheight=0.0;
      }
      pResult = Tcl_NewObj();
      z0=Location_zabs(p,pS,x0,dheight);
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(x0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(y0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(z0));
      Tcl_SetObjResult(interp, pResult);
      return TCL_OK;
    }
    pResult = Tcl_NewObj();
    for(i=2;i<objc;i+=4) {
      if( getDeckId(interp,p,objv[i],&pS) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[i+1], &x0) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(interp, objv[i+2], &y0) ) return TCL_ERROR;
      if( Tcl_GetDoubleFromObj(NULL, objv[i+3], &dheight) ) return TCL_ERROR;
      
      z0=Location_zabs(p,pS,x0,dheight);
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(x0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(y0));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(z0));
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }
  case SLICER_LOCATION_Z: {
    int i;
    double x0, y0, z0, dheight;
    Tcl_Obj *pResult;
    struct OneSlice *pS,*pAbove;
    if( objc!= 5 && objc!=6 ){
      Tcl_WrongNumArgs(interp, 2, objv, "DECK X Y ?ZOFF? ...");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pS) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &y0) ) return TCL_ERROR;
    if(objc==6) {
      if( Tcl_GetDoubleFromObj(NULL, objv[5], &dheight) ) {
        dheight=0.0;
      }
    } else {
      dheight=0.0;
    }
    z0=Location_zabs(p,pS,x0,dheight);
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(z0));
    return TCL_OK;
  }
  /*
  ** tclmethod:  SLICER actualcoords CANVAS-COORD-LIST  ?ABOVE-DECK?
  ** title:   Convert from canvas to actual coordinate space
  **
  ** The ABOVE-DECK parameter, if present, is the height above the deck.
  */
  case SLICER_ACTUALCOORDS: {
    int i, n;
    double aboveDeck;
    Tcl_Obj *pResult;
    if( objc!=3 && objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "CANVAS-COORD-LIST ?ABOVE-DECK?");
      return TCL_ERROR;
    }
    if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;
    if( n%2!=0 ){
      Tcl_AppendResult(interp, "coordinate list must contain a multiple "
         "of 2 values", 0);
      return TCL_ERROR;
    }
    if( objc==4 ){
      if( Tcl_GetDoubleFromObj(interp, objv[3], &aboveDeck) ) return TCL_ERROR;
    }else{
      aboveDeck = 0.0;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    pResult = Tcl_NewObj();
    for(i=0; i<n-1; i+=2){
      double ax, ay;
      double cx, cy;
      Tcl_Obj *pObj;
      int j;
      struct OneSlice *pS;
      Tcl_ListObjIndex(0, objv[2], i, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &cx) ) break;
      Tcl_ListObjIndex(0, objv[2], i+1, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &cy) ) break;
      for(j=0; j<p->nSlice-1 && p->a[j].upperbound>cy; j++){}
      pS = &p->a[j];

      /* Original Formula 
      ** ax = cx*p->rZoom - pS->rXShift;
      ** ay = pS->mxY + (pS->mnY-pS->mxY)*(cy-pS->top)/(pS->btm-pS->top);
      */
      ax=xCanvasToActual(p,pS,cx);
      ay=yCanvasToActual(p,pS,cy);
      
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pS->did));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(round(ax)));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(round(ay)));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(round(aboveDeck)));
    }
    if( i<n-1 ){
      Tcl_DecrRefCount(pResult);
      return TCL_ERROR;
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }
 
  /*
  ** tclmethod:  SLICER canvascoords ACTUAL-COORD-LIST ?xvar? ?yvar?
  ** title:   Convert from actual (in DID X Y) to canvas coordinate space
  **
  ** The coordinates are considered to form a line.  If the actual
  ** coordinates make a transition from one slice to another, an extra
  ** set of canvas coordinates may be inserted to show the point of
  ** this transition.
  **
  */
  case SLICER_CANVASCOORDS: {
    int i, n;
    Tcl_Obj *pResult=NULL;
    double ax, ay;
    double cx, cy;
    Tcl_Obj *pObj;
    struct OneSlice *pS;
    int error=0;

    if( objc!=3 && objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "ACTUAL-COORD-LIST ?xvar yvar?");
      return TCL_ERROR;
    }
    if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;
    if(n==3) {
      /* Goose a 3 length to a 4 length */
      n=4;
    }
    if(n<4 || n%4!=0) {
      Tcl_AppendResult(interp, "coordinate list must contain a multiple "
           "of 4 values", 0);
      return TCL_ERROR;
    }

    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if(objc != 5) {
      pResult = Tcl_NewObj();
    }

    for(i=0; i<n-3; i+=4){
      Tcl_ListObjIndex(0, objv[2], i+1, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &ax) ) {
        error=1;
        break;
      }
      Tcl_ListObjIndex(0, objv[2], i+2, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &ay) ) {
        error=1;
        break;
      }
      Tcl_ListObjIndex(0, objv[2], i, &pObj);
      if( getDeckId(interp, p, pObj, &pS) ) {
        error=1;
        break;
      }
      
      /* Old direct formula
      **
      ** cx = (ax+pS->rXShift)/p->rZoom;
      ** cy = pS->top + (pS->btm-pS->top)*(ay-pS->mxY)/(pS->mnY-pS->mxY);
      */
      cx = xActualToCanvas(p,pS,ax);       
      cy = yActualToCanvas(p,pS,ay);

      if(objc==5) {
        Tcl_ObjSetVar2(interp,objv[3],NULL,Tcl_NewDoubleObj(cx),0);
        Tcl_ObjSetVar2(interp,objv[4],NULL,Tcl_NewDoubleObj(cy),0);
        return TCL_OK;
      } else {    
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cx));
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(cy));
      }
    }
    if(error) {
      if(pResult) {
        Tcl_DecrRefCount(pResult);
      }
      return TCL_ERROR;
    }
    if( i<n-3 ){
        Tcl_AppendResult(interp, "Did not reach the end", 0);
      Tcl_DecrRefCount(pResult);
      return TCL_ERROR;
    }

    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  SLICER create DID NAME Z MIN-Y MAX-Y
  ** title:   Create a new slice
  */
  case SLICER_CREATE: {
    double z, mnY, mxY;
    const char *zNameOrig;
    char *zName;
    int did;
    int i, nName;
    if( objc!=7 ){
      Tcl_WrongNumArgs(interp, 2, objv, "DID NAME Z MIN-Y MAX-Y");
      return TCL_ERROR;
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &did) ) return TCL_ERROR;
    zNameOrig = Tcl_GetStringFromObj(objv[3], &nName);
    if( Tcl_GetDoubleFromObj(interp, objv[4], &z) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[5], &mnY) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[6], &mxY) ) return TCL_ERROR;
    if( mnY>=mxY ){
      Tcl_AppendResult(interp, "MIN-Y must be less than MAX-Y", 0);
      return TCL_ERROR;
    }
    zName = Odie_Alloc( nName+1 );
    if( zName==0 ) return TCL_ERROR;
    memcpy(zName, zNameOrig, nName+1);
    for(i=0; i<p->nSlice; i++){
      if( p->a[i].did==did ){
        Tcl_AppendResult(interp, "Deckid for slice ", zName, " is the "
           "same as existing slice ", p->a[i].zName, 0);
        return TCL_ERROR;
      }
      if( p->a[i].z==z ){
        Tcl_AppendResult(interp, "Z coordinate for slice ", zName, " is the "
           "same as existing slice ", p->a[i].zName, 0);
        return TCL_ERROR;
      }
    }
    p->nSlice++;
    p->a = Odie_Realloc((char *)p->a, sizeof(p->a[0])*p->nSlice);
    if( p->a==0 ){
      p->nSlice = 0;
      return TCL_ERROR;
    }
    for(i=p->nSlice-1; i>0 && p->a[i-1].z>z; i--){
      p->a[i] = p->a[i-1];
      p->a[i].idx = i;
    }
    p->a[i].did = did;
    p->a[i].idx = i;
    p->a[i].zName = zName;
    p->a[i].nXZ = 0;
    p->a[i].xz = 0;
    p->a[i].z = z;
    p->a[i].mnY = mnY;
    p->a[i].mxY = mxY;
    p->a[i].rXShift = p->rXOffset*z;
    computeTopAndBottom(p);
    break;
  }

  /* tclmethod:  SLICER deck X Y Z
  ** title:   Return the name of the deck at actual coordinates X,Y,Z
  **
  ** See also: find
  */
  case SLICER_DECK: {
    double x0, y0, z0;
    Tcl_Obj *pResult;
    struct OneSlice *pS;
    if( objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y Z");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &z0) ) return TCL_ERROR;
    pResult = Tcl_NewObj();
    pS = deckAt(p, x0, z0);
    pResult = ODIE_CONSTANT_STRING(pS->zName);
    Tcl_SetObjResult(interp, pResult);
    break;
  }


  /* tclmethod:  SLICER did X Y Z
  ** title:   Return the id of the deck at actual coordinates X,Y,Z
  **
  ** See also: find
  */
  case SLICER_DECKID: 
  case SLICER_DID: {
    double x0, y0, z0;
    Tcl_Obj *pResult;
    struct OneSlice *pS;
    if( objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y Z");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &z0) ) return TCL_ERROR;
    pResult = Tcl_NewObj();
    pS = deckAt(p, x0, z0);
    pResult = Tcl_NewIntObj(pS->did);
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  SLICER delete NAME
  ** title:   Remove a slice from the slicer
  */
  case SLICER_DELETE: {
    int i, j;
    const char *zName;
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME");
      return TCL_ERROR;
    }
    zName = Tcl_GetStringFromObj(objv[2], 0);
    for(i=0; i<p->nSlice && strcmp(p->a[i].zName,zName)!=0; i++){}
    if( i<p->nSlice ){
      Odie_Free((char *)p->a[i].zName);
      p->nSlice--;
      for(j=i; j<p->nSlice; j++){
        p->a[j] = p->a[j+1];
      }
      computeTopAndBottom(p);
    }
    break;
  }

  /*
  ** tclmethod:  SLICER destroy
  ** title:   Destroy this slicer
  */
  case SLICER_DESTROY: {
    Tcl_DeleteCommand(interp,Tcl_GetString(objv[0]));
    break;
  }

  case SLICER_DRAWLINE: {
    int result;
    result=slicer_method_drawline(pArg,interp,objc-1,objv+1);
    return result;
  }
  
  /*
  ** tclmethod:  SLICER objinfo obj
  ** title:   Return TK to draw a line on a canvas
  */
  case SLICER_OBJINFO: {
    Tcl_Obj *tmp;
    tmp = objv[2];
    printf("INFO: %s Ref: %d Type: %p \n", Tcl_GetStringFromObj(tmp, NULL),
        tmp->refCount, tmp->typePtr);
    fflush (stdout);
    break;
  }
  /*
  ** tclmethod:  SLICER makedlist CANVAS ACTUAL-PATH TAG TRANSDECK-TAG P-TAG ?COLOR?
  ** title:   Return TK to draw a line on a canvas
  */
  case SLICER_MAKEDLIST: {
    int i, n;
    double *aCoord;             /* Actual coordinates of line to draw */
    struct OneSlice **apDeck;   /* Array of all decks */
    Tcl_Obj *pVTag;             /* The "sNNN" tag added to all line segments */
    const char *zXTag;          /* Trans-deck tag (dashed lines) */
    Tcl_Obj *rtnList ;          /* List object to return */
    Tcl_Obj *coordList ;        /* List of coords to return */
    Tcl_Obj *configList ;       /* configuration string to return */
    Tcl_Obj *tmpObj;            /* Cheater for filling lists */

    if( objc!=5 && objc!=6 ){
      Tcl_WrongNumArgs(interp, 2, objv, 
          "ACTUAL-PATH TAG TRANSDECK-TAG ?COLOR?");
      return TCL_ERROR;
    }
// slicer makedlist  Path TagList Transdeck-TagList ?COLOR?
//  0     1          2    3       4                  5

    zXTag = Tcl_GetStringFromObj(objv[4], 0);

    if( Tcl_ListObjLength(interp, objv[2], &n) ) return TCL_ERROR;

    if( n%3!=0 || n<6 ){
      Tcl_AppendResult(interp, "coordinate list must contain a multiple "
         "of 3 values with a minimum of 6", 0);
      return TCL_ERROR;
    }

    rtnList = Tcl_NewListObj(0, NULL);

    aCoord = (double *)Odie_Alloc( sizeof(aCoord[0])*n + sizeof(apDeck[0])*(n/3) );
    if( aCoord==0 ){
      return TCL_ERROR;
    }
    
    // Move coords from Tcl objv[3] int "C" aCoord array

    apDeck = (struct OneSlice **)&aCoord[n];
    for(i=0; i<n; i++){
      Tcl_Obj *pObj;
      Tcl_ListObjIndex(0, objv[2], i, &pObj);
      if( Tcl_GetDoubleFromObj(interp, pObj, &aCoord[i]) ){
        Odie_Free((char *)aCoord);
        return TCL_ERROR;
      }
    }
    n /= 3;

    for(i=0; i<n; i++){
      double z = aCoord[i*3+2];
      double x = aCoord[i*3];
      apDeck[i] = deckAt(p, x, z);
    }

    for(i=1; i<n; i++){
      char zBuf[30];
      double x0, y0, x1, y1, z0, z1, Y0, Y1;

      coordList = Tcl_NewListObj(0, NULL);
      configList = Tcl_NewListObj(0, NULL);

      sprintf(zBuf, "s%d", i);
      pVTag = Tcl_NewStringObj(zBuf, -1);
      if( apDeck[i]!=apDeck[i-1] ) {

        /* Old direct formula
        **
        ** x0=(aCoord[i*3-3]+apDeck[i-1]->rXShift)/p->rZoom;
        ** x1 = (aCoord[i*3]+apDeck[i]->rXShift)/p->rZoom;
        */
        x0 = xActualToCanvas(p,apDeck[i-1],aCoord[i*3-3]);
        x1 = xActualToCanvas(p,apDeck[i],aCoord[i*3]);
        y0 = yActualToCanvas(p,apDeck[i-1], Y0=aCoord[i*3-2]);
        y1 = yActualToCanvas(p,apDeck[i], Y1=aCoord[i*3+1]);
        z0 = aCoord[i*3-1];
        z1 = aCoord[i*3+2];
        if( zXTag[0]!=0 ){
          Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(x0));
          Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(y0));
          Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(x1));
          Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(y1));

          Tcl_ListObjAppendElement(interp, configList,
	      ODIE_CONSTANT_STRING("-tags"));

          tmpObj = Tcl_DuplicateObj(objv[4]);
          Tcl_ListObjAppendElement(interp, tmpObj, pVTag);
          Tcl_ListObjAppendElement(interp, configList, tmpObj);

        }
      } else {
        /* Old direct formula
        **
        ** x0 = (aCoord[i*3-3]+apDeck[i]->rXShift)/p->rZoom;
        ** x1 = (aCoord[i*3]+apDeck[i]->rXShift)/p->rZoom;
        */
        x0 = xActualToCanvas(p,apDeck[i],aCoord[i*3-3]);
        x1 = xActualToCanvas(p,apDeck[i],aCoord[i*3]);
        
        y0 = yActualToCanvas(p,apDeck[i], aCoord[i*3-2]);
        y1 = yActualToCanvas(p,apDeck[i], aCoord[i*3+1]);

        Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(x0));
        Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(y0));
        Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(x1));
        Tcl_ListObjAppendElement(interp, coordList,Tcl_NewDoubleObj(y1));


          Tcl_ListObjAppendElement(interp, configList,\
	      ODIE_CONSTANT_STRING("-tags"));

          tmpObj = Tcl_DuplicateObj(objv[3]);
          Tcl_ListObjAppendElement(interp, tmpObj, pVTag);

          Tcl_ListObjAppendElement(interp, configList, tmpObj);

      }
      if( objc>=6 ) {
         Tcl_ListObjAppendElement(interp, configList,\
           ODIE_CONSTANT_STRING("-fill"));
         Tcl_ListObjAppendElement(interp, configList,\
	   Tcl_DuplicateObj(objv[5]));
      }


    Tcl_ListObjAppendElement(interp, rtnList, coordList);
    Tcl_ListObjAppendElement(interp, rtnList, configList);
    }
    Odie_Free((char *)aCoord);

    Tcl_SetObjResult(interp, rtnList);

    break;
  }

  /* tclmethod:  SLICER find X Y
  ** title:   Return the name of the deck at canvas coordinates X,Y
  **
  ** The "deck" command works similarly except that it uses actual
  ** coordinates as inputs.
  */
  case SLICER_FIND: {
    int i;
    double x0, y0;
    Tcl_Obj *pResult;
    struct OneSlice *pS;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    pResult = Tcl_NewObj();
    for(i=0; i<p->nSlice-1 && p->a[i].upperbound>y0; i++){}
    pS = &p->a[i];
    pResult = ODIE_CONSTANT_STRING(pS->zName);
    Tcl_SetObjResult(interp, pResult);
    break;
  }

 /* tclmethod:  SLICER finddid X Y
  ** title:   Return the name of the deck at canvas coordinates X,Y
  **
  ** The "deck" command works similarly except that it uses actual
  ** coordinates as inputs.
  */
  case SLICER_FINDDID: {
    int i;
    double x0, y0;
    Tcl_Obj *pResult;
    struct OneSlice *pS;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    pResult = Tcl_NewObj();
    for(i=0; i<p->nSlice-1 && p->a[i].upperbound>y0; i++){}
    pS = &p->a[i];
    pResult = Tcl_NewIntObj(pS->did);
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /* tclmethod:  SLICER height NAME X
  ** title:   Return the height of slice NAME at actual position X
  **            Assuming the deck were flat
  **
  ** See also "headroom".
  */
  case  SLICER_FLATHEIGHT: {
    struct OneSlice *ppS;
    double x0;
    
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME X");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;
    if( getDeckId(interp,p,objv[2],&ppS) ) return TCL_ERROR;

    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(ppS->z));
    break;
  }

  /* tclmethod:  SLICER headroom NAME X
  ** title:   Return the headroom (Z from deck) of slice NAME at actual position X
  **
  ** See also "height"
  */
  case SLICER_HEADROOM: {
    double x0;
    struct OneSlice *pDeck, *pAbove;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME X");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;

    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    pAbove = deckAbove(p, pDeck);
    if( pAbove ){
      Tcl_SetObjResult(interp,Tcl_NewDoubleObj(deckHeight(pAbove,x0)-deckHeight(pDeck,x0)));
    }else{
      Tcl_SetObjResult(interp,Tcl_NewDoubleObj(p->upper_height));
    }
    break;
  }

  /* tclmethod:  SLICER ceiling NAME X ?DEFAULT?
  ** title:   Return the Z (absolute) of the ceiling of slice NAME at actual position X
  **
  ** See also "height"
  */
  case SLICER_CEILING: {
    double x0;
    struct OneSlice *pDeck, *pAbove;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME X");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;

    if( getDeckId(interp,p,objv[2],&pDeck) ) return TCL_ERROR;
    pAbove = deckAbove(p, pDeck);
    if( pAbove ){
      Tcl_SetObjResult(interp,Tcl_NewDoubleObj(deckHeight(pAbove,x0)));
    }else{
      Tcl_SetObjResult(interp,Tcl_NewDoubleObj(deckHeight(pDeck,x0)+p->upper_height));
    }
    break;
  }
  
  /* tclmethod:  SLICER height NAME X ?zoff?
  ** title:   Return the height (absolute) of slice NAME at actual position X
  ** description:
  **   NOTE if the top deck is asked for, and zoff is negative the
  **   system will assume a 2000mm ceiling
  ** See also "headroom".
  */
  case  SLICER_HEIGHT: {
    struct OneSlice *ppS,*pAbove;
    double x0,zoff=0,zresult=0.0;
    if( objc!=4 && objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME X ?ZOFF?");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;
    if( getDeckId(interp,p,objv[2],&ppS) ) return TCL_ERROR;
    if( objc==5 ) {
      if( Tcl_GetDoubleFromObj(NULL, objv[4], &zoff) ) {
        zoff=0.0;
      }
    }
    if(zoff < 0.0) {
      pAbove = deckAbove(p, ppS);
      if(pAbove) {
        zresult=deckHeight(pAbove,x0)+zoff;
      } else {
        zresult=deckHeight(ppS,x0)+2000+zoff;
      }
    } else {
      zresult=deckHeight(ppS,x0)+zoff;
    }
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(zresult));
    break;
  }

  /* tclmethod:  SLICER deckheight NAME X ?zoff?
  ** title:   Return the height (mm from floor) of slice NAME at actual position X
  ** description:
  **   NOTE if the top deck is asked for, and zoff is negative the
  **   system will assume a 2000mm ceiling
  ** See also "headroom".
  */
  case  SLICER_DECKHEIGHT: {
    int i;
    double x0, z0, dheight;
    Tcl_Obj *pResult;
    struct OneSlice *pS,*pAbove;
    if( objc!= 4 && objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "DECK X ?ZOFF?");
      return TCL_ERROR;
    }
    if( p->nSlice<=0 ){
      Tcl_AppendResult(interp, "no slices defined", 0);
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&pS) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x0) ) return TCL_ERROR;
    if(objc==5) {
      if( Tcl_GetDoubleFromObj(NULL, objv[4], &dheight) ) {
        dheight=0.0;
      }
    } else {
      dheight=0.0;
    }
    z0=Location_zdeck(p,pS,x0,dheight);
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(z0));
    return TCL_OK;
  }

  /*
  ** tclmethod:  SLICER info NAME
  ** title:   Return information about a particular slice
  */
  case SLICER_INFO: {
    struct OneSlice *ppS;
    
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME");
      return TCL_ERROR;
    }
    if( getDeckId(interp,p,objv[2],&ppS) ) {
      return TCL_ERROR;
    } else {
      Tcl_Obj *pResult = Tcl_NewObj();
      
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("name"));
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING(ppS->zName));  
 
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("did"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewIntObj(ppS->did));  
 
    
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("z"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->z));

      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("miny"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->mnY));
      
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("maxy"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->mxY));
      
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("top"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->top));
      
      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("bottom"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->btm));

      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("above"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->above));

      Tcl_ListObjAppendElement(interp, pResult, ODIE_CONSTANT_STRING("below"));
      Tcl_ListObjAppendElement(interp, pResult, Tcl_NewDoubleObj(ppS->below));
      
      Tcl_SetObjResult(interp, pResult);
    }
    break;
  }

  /* 
  ** tclmethod:  SLICER list
  ** title:   List the names of all defined slices in the slicer
  */
  case SLICER_LIST: {
    Tcl_Obj *pResult = Tcl_NewObj();
    int i;
    for(i=0; i<p->nSlice; i++){
      Tcl_ListObjAppendElement(0, pResult, ODIE_CONSTANT_STRING(p->a[i].zName));
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }


  /*
  ** tclmethod:  SLICER profile DECK ?X Z ...?
  ** title:   Create an inboard profile for a deck
  */
  case SLICER_PROFILE: {
    int i, j, min;
    struct OneSlice *pS;
    if( objc<3 || (objc>3 && objc<7) || (objc & 1)==0 ){
      Tcl_WrongNumArgs(interp, 2, objv, "NAME ?X Z X Z...?");
      return TCL_ERROR;
    }
    
    if(getDeckId(interp,p,objv[2],&pS)) {
      return TCL_ERROR;
    }
    if( objc==3 ){
      Tcl_Obj *pResult = Tcl_NewObj();
      for(i=0; i<pS->nXZ; i++){
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(pS->xz[i]));
      }
      Tcl_SetObjResult(interp, pResult);
    }else{
      double *xz;
      if( pS->xz ) Odie_Free((char *)pS->xz);
      pS->nXZ = objc - 3;
      xz = pS->xz = (double *)Odie_Alloc( sizeof(pS->xz[0])*pS->nXZ );
      if( xz==0 ){
        pS->nXZ = 0;
        break;
      }
      for(i=0; i<pS->nXZ; i++){
        if( Tcl_GetDoubleFromObj(interp, objv[i+3], &xz[i]) ){
          Odie_Free((char *)xz);
          pS->nXZ = 0;
          pS->xz = 0;
          return TCL_ERROR;
        }
      }

      /* Put the profile in increasing X order.  An N**2 sort is used because
      ** it is convenient and the list will usually be short. */
      for(i=0; i<pS->nXZ-2; i+=2){
        for(min=i, j=i+2; j<pS->nXZ; j+=2){
          if( xz[j]<xz[min] ) min = j;
        }
        if( min>i ){
          double t = xz[min];
          xz[min] = xz[i];
          xz[i] = t;
          t = xz[min+1];
          xz[min+1] = xz[i+1];
          xz[i+1] = t;
        }
      }

      /* Remove duplidate X coordinates */
      for(i=j=0; i<pS->nXZ; i+=2){
        if( i<pS->nXZ-2 && xz[i+2]==xz[i] ){
          /* Ignore the duplicate */
        }else{
          xz[j++] = xz[i];
          xz[j++] = xz[i+1];
        }
      }
      pS->nXZ = j;
      if( j<4 ){
        pS->nXZ = 0;
        Odie_Free((char *)pS->xz);
        pS->xz = 0;
      }
    }        
    break;
  }

  /*
  ** tclmethod:  SLICER xoffset ?AMT?
  ** title:   Change the X-Offset as a function of deck height
  */
  case SLICER_XOFFSET: {
    double rXOffset;
    int i;
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==2 ){
      Tcl_SetObjResult(interp, Tcl_NewDoubleObj(p->rXOffset));
      return TCL_OK;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &rXOffset) ) return TCL_ERROR;
    for(i=0; i<p->nSlice; i++){
      p->a[i].rXShift = rXOffset*p->a[i].z;
    }
    p->rXOffset = rXOffset;
    break;
  }

  /*
  ** tclmethod:  SLICER zoom ?ZOOM?
  ** title:   Query or change the zoom factor.
  */
  case SLICER_ZOOM: {
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==3 ){
      double r;
      if( Tcl_GetDoubleFromObj(interp, objv[2], &r) ) return TCL_ERROR;
      p->rZoom = r;
    }
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(p->rZoom));
    computeTopAndBottom(p);
    break;
  }

  /*
  ** tclmethod:  SLICER upper_height MM
  ** title:   Set the headroom on the top level of the slicer (defaults to 2000)
  */
  case SLICER_UPPER_HEIGHT: {
    if( objc != 2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "Z");
      return TCL_ERROR;
    }
    if( objc==3 ){
      double r;
      if( Tcl_GetDoubleFromObj(interp, objv[2], &r) ) return TCL_ERROR;
      p->upper_height=r;
    }
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(p->upper_height));
    break;
  }

  /* End of the command methods.  The brackets that follow terminate the
  ** automatically generated switch.
  ****************************************************************************/
  }
  }
  return TCL_OK;
}

/*
** tclcmd: slicer SLICER
** title: creates a slicer object
** This routine runs when the "slicer" command is invoked to create a
** new slicer.
*/
int Odie_SlicerCreateProc(
  void *NotUsed,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  char *zCmd;
  Slicer *p;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "SLICER");
    return TCL_ERROR;
  }
  zCmd = Tcl_GetStringFromObj(objv[1], 0);
  p = (Slicer *)Odie_Alloc( sizeof(*p) );
  p->rZoom = 100;
  p->upper_height=2000;
  Tcl_CreateObjCommand(interp, zCmd, slicerMethodProc, p, destroySlicer);
  return TCL_OK;
}

/*** Automatically Generated Header File - Do Not Edit ***/
  const static char *SLICER_strs[] = {
    "above",              "actualcoords",      "below",
    "canvascoords",       "ceiling",           "create",
    "deck",               "deckheight",        "deckid",
    "deckid_to_name",     "deckname_to_id",    "delete",
    "destroy",            "did",               "drawline",
    "find",               "finddid",           "flatheight",
    "headroom",           "height",            "info",
    "list",               "location_to_xyz",   "location_z",
    "makedlist",          "objinfo",           "profile",
    "upper_height",       "xoffset",           "xyz_to_location",
    "zoom",               0                    
  };
  enum SLICER_enum {
    SLICER_ABOVE,         SLICER_ACTUALCOORDS, SLICER_BELOW,
    SLICER_CANVASCOORDS,  SLICER_CEILING,      SLICER_CREATE,
    SLICER_DECK,          SLICER_DECKHEIGHT,   SLICER_DECKID,
    SLICER_DECKID_TO_NAME, SLICER_DECKNAME_TO_ID,SLICER_DELETE,
    SLICER_DESTROY,       SLICER_DID,          SLICER_DRAWLINE,
    SLICER_FIND,          SLICER_FINDDID,      SLICER_FLATHEIGHT,
    SLICER_HEADROOM,      SLICER_HEIGHT,       SLICER_INFO,
    SLICER_LIST,          SLICER_LOCATION_TO_XYZ,SLICER_LOCATION_Z,
    SLICER_MAKEDLIST,     SLICER_OBJINFO,      SLICER_PROFILE,
    SLICER_UPPER_HEIGHT,  SLICER_XOFFSET,      SLICER_XYZ_TO_LOCATION,
    SLICER_ZOOM,         
  };
 int index;
  if( objc<2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "METHOD ?ARG ...?");
    return TCL_ERROR;
  }
  if( Tcl_GetIndexFromObj(interp, objv[1], SLICER_strs, "option", 0, &index)){
    return TCL_ERROR;
  }
  switch( (enum SLICER_enum)index )

/*
** This file implements a TCL object that keeps track of the walls and
** bulkheads on a single deck of a ship.
**
** This widget assumes a right-handed coordinate system if zoom is positive
** and a left-handed coordinate system is zoom is negative.  The Tk canvas
** widget uses a left-handed coordinate system all the time.  The READI
** database uses a right-handed coordinate system all the time.  This module
** can be used to translate by setting zoom to +1.0 for database I/O and
** to -$g(zoom) for canvas I/O.
**
** This module uses a purely 2-D model.  It can only handle a single
** deck at a time.  If a multi-deck model needs to be displayed then
** that multi-deck model should first be flattened into a stack of
** individual decks in the same plane using the separate "slicer" object.
**
** This file implements a single new constructor tcl command named "wallset".
** The wallset command creates a new wallset object.  Methods on this
** wallset object are used to manage the object.
**
** The details of the various methods and what they do are provided in
** header comments above the implementation of each method.
*/
const char wallset_c_version[] = "$Header: /readi/code/tobe/wallset.c,v 1.34 2007/02/22 15:04:04 drh Exp $";
#include "odieInt.h"
#include <stdarg.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <math.h>

#ifndef M_PI
# define M_PI 3.1415926535898
#endif

/*
** Remove all of the ComptBox entries from the wallset.
*/
static void clearComptBoxCache(Wallset *pWS){
  ComptBox *p = pWS->pComptBox;
  while( p ){
    ComptBox *pNext = p->pNext;
    Odie_Free((char *)p);
    p = pNext;
  }
  pWS->pComptBox = 0;
}

/*
** This routine is invoked when the TCL command that implements a
** wallset is deleted.  Free all memory associated with that
** wallset.
*/
static void destroyWallset(void *pArg){
  Wallset *p = (Wallset*)pArg;
  Link *pLink = p->pAll;
  clearComptBoxCache(p);
  while( pLink ){
    Segment *pSeg = pLink->pLinkNode;
    pLink = pSeg->pAll.pNext;
    Odie_Free((char *) pSeg );
  }
  Odie_Free((char *) p );
}

/*
** Clear the Segment.ignore flag on all segments within a wallset.
*/
static void ignoreNone(Wallset *p){
#if 0  
  Link *pLink;
  for(pLink=p->pAll; pLink; pLink=pLink->pNext){
    pLink->pSeg->ignore = 0;
  }
#endif
}

/*
** Return a pointer to the segment with the given ID.  Return NULL
** if there is no such segment.
*/
static Segment *findSegment(Wallset *p, int id){
  int h;
  Link *pLink;

  h = hashInt(id);
  for(pLink = p->hashId[h]; pLink; pLink=pLink->pNext){
    Segment *pSeg=pLink->pLinkNode;
    if( pSeg->id==id ) return pSeg;
  }
  return 0;
}

#if 0 /* NOT USED */
/*
** Make a copy of a string.
*/
static char *stringDup(const char *z){
  int len = strlen(z);
  char *zNew = Odie_Alloc( len+1 );
  if( zNew ){
    memcpy(zNew, z, len+1);
  }
  return zNew;
}
#endif

/*
** Scan all segments looking for the vertex or vertices that are nearest
** to x,y.  Return a pointer to a Segment.set that is the list of matching
** segments.  Also write the nearest point into *pX,*pY.
**
** The returned list uses the Segment.set link.
*/
static Link *nearestVertex(
  Wallset *p,               /* The wallset to be scanned */
  double x, double y,       /* Search for points near to this point */
  double *pX, double *pY    /* Write nearest vertex here */
){
  double nx, ny;
  double min = -1.0;
  Link *pList = 0;
  Link *pI;

  x = roundCoord(x);
  y = roundCoord(y);
  for(pI=p->pAll; pI; pI=pI->pNext){
    double dx, dy, dist;
    Segment *pSeg = pI->pLinkNode;
    dx = x - pSeg->from[X_IDX];
    dy = y - pSeg->from[Y_IDX];
    dist = dx*dx + dy*dy;
    if( min<0.0 || dist<=min ){
      if( min<0.0 || nx!=pSeg->from[X_IDX] || ny!=pSeg->from[Y_IDX] ){
        pList = 0;
        nx = pSeg->from[X_IDX];
        ny = pSeg->from[Y_IDX];
        min = dist;
      }
      LinkInit(pSeg->pSet, pSeg);
      LinkInsert(&pList, &pSeg->pSet);
    }
    dx = x - pSeg->to[X_IDX];
    dy = y - pSeg->to[Y_IDX];
    dist = dx*dx + dy*dy;
    if( dist<=min ){
     if( nx!=pSeg->to[X_IDX] || ny!=pSeg->to[Y_IDX] ){
        pList = 0;
        nx = pSeg->to[X_IDX];
        ny = pSeg->to[Y_IDX];
        min = dist;
      }
      LinkInit(pSeg->pSet, pSeg);
      LinkInsert(&pList, &pSeg->pSet);
    }
  }
  *pX = nx;
  *pY = ny;
  return pList;
}

/*
** Scan all segments looking for the point on a segment that is nearest
** to x,y.  Return a pointer to a Segment.set that is the list of matching
** segments.  This set might contain multiple members if the nearest point
** is actually a vertex shared by two or more segments.  Write the nearest
** point into *pX, *pY.
**
** /// Ignore any segment that has its Segment.ignore flag set. -- removed
**
** The returned list uses the Segment.set list.
*/
static Link *nearestPoint(
  Wallset *p,               /* The wallset to be scanned */
  double x, double y,       /* Search for points near to this point */
  double *pX, double *pY    /* Write nearest vertex here */
){
  double nx, ny;
  double min = -1.0;
  Link *pList = 0;
  Link *pI;

  x = roundCoord(x);
  y = roundCoord(y);
  for(pI=p->pAll; pI; pI=pI->pNext){
    double dx, dy, dist;
    Segment *pSeg;
    double acx, acy;   /* Vector from x0,y0 to x,y */
    double abx, aby;   /* Vector from x0,y0 to x1,y1 */
    double rx, ry;     /* Nearest point on x0,y0->to[X_IDX],y1 to x,y */
    double r;

    pSeg = pI->pLinkNode;
    /* if( pSeg->ignore ) continue; */
    acx = x - pSeg->from[X_IDX];
    acy = y - pSeg->from[Y_IDX];
    abx = pSeg->to[X_IDX] - pSeg->from[X_IDX];
    aby = pSeg->to[Y_IDX] - pSeg->from[Y_IDX];
    r = (acx*abx + acy*aby)/(abx*abx + aby*aby);
    if( r<=0 ){
      rx = pSeg->from[X_IDX];
      ry = pSeg->from[Y_IDX];
    }else if( r>=1 ){
      rx = pSeg->to[X_IDX];
      ry = pSeg->to[Y_IDX];
    }else{
      rx = pSeg->from[X_IDX] + abx*r;
      ry = pSeg->from[Y_IDX] + aby*r;
    }
    rx = roundCoord(rx);
    ry = roundCoord(ry);
    dx = x - rx;
    dy = y - ry;
    dist = dx*dx + dy*dy;
    if( min<0.0 || dist<=min ){
      if( min<0.0 || nx!=rx || ny!=ry ){
        pList = 0;
        nx = rx;
        ny = ry;
        min = dist;
      }
      LinkInit(pSeg->pSet, pSeg);
      LinkInsert(&pList, &pSeg->pSet);
    }
  }
  *pX = nx;
  *pY = ny;
  return pList;
}

/*
** Return TRUE if the value x is in between x1 and x2.
*/
static int between(double x, double x1, double x2){
  if( x1<x2 ){
    return x>=x1 && x<=x2;
  }else{
    return x>=x2 && x<=x1;
  }
}

/*
** Return TRUE if the given segment is on the given list
*/
static int segmentOnList(Segment *pSeg, Link *pList){
  while( pList ){
    if( pList->pLinkNode==pSeg ) return 1;
    pList = pList->pNext;
  }
  return 0;
}

/*
** Return a list of all segments which have an end at the given vertex.
** The returned list uses Segment.set
*/
static Link *segmentsAtVertex(Wallset *p, double x, double y){
  Link *pList = 0;
  Link *pI;
  int h;

  x = roundCoord(x);
  y = roundCoord(y);
  h = hashCoord(x, y);
  for(pI=p->hashFrom[h]; pI; pI=pI->pNext){
    Segment *pSeg = pI->pLinkNode;
    /* if( pSeg->ignore ) continue; */
    if( floatCompare(x, pSeg->from[X_IDX])==0 && floatCompare(y, pSeg->from[Y_IDX])==0 ){
      assert( !segmentOnList(pSeg, pList) );
      LinkInit(pSeg->pSet, pSeg);
      LinkInsert(&pList, &pSeg->pSet);
    }
  }
  for(pI=p->hashTo[h]; pI; pI=pI->pNext){
    Segment *pSeg = pI->pLinkNode;
    /* if( pSeg->ignore ) continue; */
    if( floatCompare(x, pSeg->to[X_IDX])==0 && floatCompare(y, pSeg->to[Y_IDX])==0 ){
      assert( !segmentOnList(pSeg, pList) );
      LinkInit(pSeg->pSet, pSeg);
      LinkInsert(&pList, &pSeg->pSet);
    }
  }
  return pList;
}

/*
** The point xV,yV is a vertex in the wallset.  This routine locates
** a segment connected to that vertex which is the first segment in
** a clockwise direction from xR,yR->xV,yV.  A pointer to the segment
** is written into *ppSeg.  If the output segment moves backwards
** (in other words if x1,y1 of the segment is connected at xV,yV)
** then *pfBack is true.
**
** If a suitable segment is found, 0 is returned.  Non-zero is returned
** if no suitable segment could be found.
**
** This routine uses the Segment.set list internally.
*/
static int nextCwSegment(
  Wallset *p,              /* The wallset */
  double xR, double yR,    /* Remote end of input segment */
  double xV, double yV,    /* Vertex (near end of input segment) */
  Segment **ppSeg,         /* OUT: First segment clockwise from xR,yR->xV,yV */
  int *pfBack              /* OUT: True if output segment goes backwards */
){
  Link *pList, *pI;
  double rRef, rBest;
  int i, nSeg, iBest;
  Segment *pSeg;
  struct {
    Segment *pSeg;
    int isBack;
    double rAngle;
  } *aSeg, aSegStatic[20];

  /* Find all segments at xV,yV */
  pList = segmentsAtVertex(p, xV, yV);
  for(pI=pList, nSeg=0; pI; nSeg++, pI=pI->pNext){}
  if( nSeg==0 ) return 1;
  if( nSeg<=sizeof(aSegStatic)/sizeof(aSegStatic[0]) ){
    aSeg = aSegStatic;
  }else{
    aSeg = (void *)Odie_Alloc( nSeg*sizeof(*aSeg) );
  }
  for(pI=pList, i=0; pI; i++, pI=pI->pNext){ 
    aSeg[i].pSeg = pSeg = pI->pLinkNode;
    aSeg[i].isBack = floatCompare(xV, pSeg->to[X_IDX])==0 
                          && floatCompare(yV, pSeg->to[Y_IDX])==0;
  }
  
  /* Find the reference angle */
  rRef = atan2(yR-yV, xR-xV)*180.0/M_PI;

  /* Find angles on all segments */
  for(i=0; i<nSeg; i++){
    pSeg = aSeg[i].pSeg;
    if( aSeg[i].isBack ){
      aSeg[i].rAngle = atan2(pSeg->from[Y_IDX]-pSeg->to[Y_IDX], pSeg->from[X_IDX]-pSeg->to[X_IDX])*180.0/M_PI;
    }else{
      aSeg[i].rAngle = atan2(pSeg->to[Y_IDX]-pSeg->from[Y_IDX], pSeg->to[X_IDX]-pSeg->from[X_IDX])*180.0/M_PI;
    }
  }

  /* Subtract 360 to any segment angle that is less than the reference angle */
  for(i=0; i<nSeg; i++){
    if( aSeg[i].rAngle<rRef ) aSeg[i].rAngle += 360;
  }

  /* Choose the segment with the largest angle */
  rBest = aSeg[0].rAngle;
  iBest = 0;
  for(i=1; i<nSeg; i++){
    if( aSeg[i].rAngle>rBest ){
      iBest = i;
      rBest = aSeg[i].rAngle;
    }
  }
  *ppSeg = aSeg[iBest].pSeg;
  *pfBack = aSeg[iBest].isBack;
  if( aSeg!=aSegStatic ){
    Odie_Free((char *) aSeg );
  }

  return 0;
}

/*
** Consider a line beginning at x0,y0 then going from x1,y1 to x2,y2.
** x1,y1 is an elbow in the line.  This routine returns -1 if the
** elbow bends to the right, and +1 if it bends to the left.  zero is
** returned if the elbow does not bend at all.
*/
static int bendDirection(
  double x0, double y0,
  double x1, double y1,
  double x2, double y2
){
  /* Algorithm:  Rotate x0,y0->to[X_IDX],y1 90 degrees counter-clockwise.  Take
  ** the dot product with x1,y1->x2,y2.  The dot produce will be the product
  ** of two (non-negative) magnitudes and the cosine of the angle.  So if
  ** the dot product is positive, the bend is to the left, or to the right if
  ** the dot product is negative.
  */
  double r = (y0-y1)*(x2-x1) + (x1-x0)*(y2-y1);
  return r<0.0 ? +1 : (r>0.0 ? -1 : 0);
}

/*
** Given an interior point xI,yI, this routine finds a segment on the
** boundary that contains the interior point.  That segment is returned
** in *ppSeg.  *pfLeft is set to true if the interior point is to the left
** of the segment and false if it is to the right.
**
** Zero is returned on success.  Non-zero is returned if no suitable
** boundary could be located.  Non-zero might be returned, for example,
** if xI,yI is positioned directly on top of a wall or if there are no
** walls in the wallset.
**
** // Any segment marked with Segment.ignore is ignored for purposes of
** // this routine.  -- removed
**
** This routine uses the Segment.set list internally.
*/
static int firstBoundarySegment(
  Wallset *p,              /* The wallset */
  double xI, double yI,    /* An interior point */
  Segment **ppSeg,         /* OUT: A segment on the boundary containing xI,yI */
  int *pfLeft              /* OUT: True if xI,yI is to the left side *ppSeg */
){
  Link *pList;
  double xN, yN;

  /* Find nearest point, xN,yN */
  pList = nearestPoint(p, xI, yI, &xN, &yN);
  if( pList==0 ) return 1;
  if( pList->pNext ){
    /* xN,yN is a vertex...
    ** Locate the first segment clockwise from xI,yI->xN,yN and return
    */
    return nextCwSegment(p, xI, yI, xN, yN, ppSeg, pfLeft);
  }else{
    /* xN,yN is a point on single line segment...
    */
    Segment *pSeg;
    pSeg = *ppSeg = pList->pLinkNode;
    *pfLeft = bendDirection(pSeg->from[X_IDX], pSeg->from[Y_IDX], xN, yN, xI, yI)>0;
  }
  return 0;
}

/*
** Fill the given Boundary array with a list of segments (with
** Segment.ignore set to false) that form a closed circuit.  The
** first entry in aBound[] has already been filled in by the
** calling function and is used to seed the search.
**
** At most nBound slots in aBound[] will be used.  The return value
** is the number of slots in aBound[] that would have been used if those
** slots had been available.  A return of 0 indicates that no boundary
** is available.
**
** If the checkIsPrimary flag is true and the aBound[0] entry is not
** the primary segment for the compartment, then the aBound[] is not
** completely filled in and the routine returns 0;
*/
static int completeBoundary(
  Wallset *p,             /* The wallset */
  int checkIsPrimary,     /* Abort if aBound[0] is not the primary segment */
  int nBound,             /* Number of slots available in aBound[] */
  Boundary *aBound        /* IN-OUT: Write results into aBound[1...] */
){
  int cnt = 1;
  Segment *pSeg, *pS;
  int isLeft;
  int isBack;
  double xR, yR, xV, yV;

  pS = pSeg = aBound[0].pSeg;
  isLeft = aBound[0].backwards;
  if( !isLeft ){
    xR = pSeg->from[X_IDX];
    yR = pSeg->from[Y_IDX];
    xV = pSeg->to[X_IDX];
    yV = pSeg->to[Y_IDX];
  }else{
    xV = pSeg->from[X_IDX];
    yV = pSeg->from[Y_IDX];
    xR = pSeg->to[X_IDX];
    yR = pSeg->to[Y_IDX];
  }
  while( nextCwSegment(p,xR,yR,xV,yV,&pS,&isBack)==0 &&
              (isBack!=isLeft || pS!=pSeg) ){
    if( checkIsPrimary ){
      if( pS->id<pSeg->id ) return 0;
      if( pS->id==pSeg->id && !isLeft ) return 0;
    }
    if( isBack ){
      xV = pS->from[X_IDX];
      yV = pS->from[Y_IDX];
      xR = pS->to[X_IDX];
      yR = pS->to[Y_IDX];
    }else{
      xR = pS->from[X_IDX];
      yR = pS->from[Y_IDX];
      xV = pS->to[X_IDX];
      yV = pS->to[Y_IDX];
    }
    if( nBound>cnt ){ 
      aBound[cnt].pSeg = pS;
      aBound[cnt].backwards = isBack;
    }
    cnt++;
    if( cnt>1000 /* 00 */ ) return -cnt;   /* Avoid an infinite loop */
  }
  return cnt;
}

/*
** Compute the "spin" on a boundary.  A positive value means the
** circulation is to counter-clockwise and a negative value means the 
** circulation is clockwise.  For boundaries, a positive
** value means the region is internal and a negative value means
** the region is external.
*/
static double spin(Boundary *aBound, int nBound){
  double sum = 0;
  int i;
  for(i=0; i<nBound; i++){
    double x0, y0, x1, y1;
    double dx, dy;
    Segment *pSeg = aBound->pSeg;
    if( aBound->backwards ){
      x0 = pSeg->to[X_IDX];
      y0 = pSeg->to[Y_IDX];
      x1 = pSeg->from[X_IDX];
      y1 = pSeg->from[Y_IDX];
    }else{
      x0 = pSeg->from[X_IDX];
      y0 = pSeg->from[Y_IDX];
      x1 = pSeg->to[X_IDX];
      y1 = pSeg->to[Y_IDX];
    }
    aBound++;
    dx = x1-x0;
    dy = y1-y0;
    sum += x0*dy - y0*dx;
  }
  return sum;
}

/*
** The input is two linked lists of ComptBox structures where each
** list is sorted by increasing area.  Combine these two lists into
** a single sorted linked list.
*/
static ComptBox *mergeComptBox(ComptBox *p1, ComptBox *p2){
  ComptBox head;
  ComptBox *pTail = &head;
  ComptBox *p;
  while( p1 && p2 ){
    if( p1->area<=p2->area ){
      p = p1->pNext;
      pTail->pNext = p1;
      pTail = p1;
      p1 = p;
    }else{
      p = p2->pNext;
      pTail->pNext = p2;
      pTail = p2;
      p2 = p;
    }
  }
  if( p1 ){
    pTail->pNext = p1;
  }else{
    pTail->pNext = p2;
  }
  return head.pNext;
}

/*
** Construct the ComptBox cache.  For each compartment (where a compartment
** is a closed circuit of Segments) make an entry on the Wallset.pComptBox
** list.
**
** If the ComptBox cache already exists, this routine is a no-op.
*/
static void buildComptBoxCache(Wallset *p){
  Link *pI;
  int i;
  ComptBox *aSort[30];

  /* Return immediately if the cache already exists */
  if( p->pComptBox ) return;

  /* Compute a linked list of all compartment boxes */
  for(pI=p->pAll; pI; pI=pI->pNext){
    int i, j, n;
    Boundary aBound[1000];

    aBound[0].pSeg = pI->pLinkNode;
    for(j=0; j<2; j++){
      aBound[0].backwards = j;
      n = completeBoundary(p, 1, sizeof(aBound)/sizeof(aBound[0]), aBound);
      if( n>0 && spin(aBound,n)>0.0 ){
        double dx, dy;
        Segment *pSeg = pI->pLinkNode;
        ComptBox *pNew = (ComptBox *)Odie_Alloc( sizeof(*pNew) );
        pNew->pNext = p->pComptBox;
        pNew->bbox.l = pNew->bbox.r = pSeg->from[X_IDX];
        pNew->bbox.t = pNew->bbox.b = pSeg->from[Y_IDX];
        pNew->prim = aBound[0];
        for(i=1; i<n; i++){
          Segment *pSeg = aBound[i].pSeg;
          if( pSeg->from[X_IDX]<pNew->bbox.l ) pNew->bbox.l = pSeg->from[X_IDX];
          if( pSeg->from[X_IDX]>pNew->bbox.r ) pNew->bbox.r = pSeg->from[X_IDX];
          if( pSeg->from[Y_IDX]<pNew->bbox.b ) pNew->bbox.b = pSeg->from[Y_IDX];
          if( pSeg->from[Y_IDX]>pNew->bbox.t ) pNew->bbox.t = pSeg->from[Y_IDX];
          if( pSeg->to[X_IDX]<pNew->bbox.l ) pNew->bbox.l = pSeg->to[X_IDX];
          if( pSeg->to[X_IDX]>pNew->bbox.r ) pNew->bbox.r = pSeg->to[X_IDX];
          if( pSeg->to[Y_IDX]<pNew->bbox.b ) pNew->bbox.b = pSeg->to[Y_IDX];
          if( pSeg->to[Y_IDX]>pNew->bbox.t ) pNew->bbox.t = pSeg->to[Y_IDX];
        }
        dx = pNew->bbox.r - pNew->bbox.l;
        dy = pNew->bbox.t - pNew->bbox.b;
        pNew->area = sqrt(dx*dx+dy*dy);
        p->pComptBox = pNew;
      }
    }
  }

  /* Sort the list into order of increasing area */
  for(i=0; i<sizeof(aSort)/sizeof(aSort[0]); i++) aSort[i] = 0;
  while( p->pComptBox ){
    ComptBox *pBox = p->pComptBox;
    p->pComptBox = pBox->pNext;
    pBox->pNext = 0;
    for(i=0; i<sizeof(aSort)/sizeof(aSort[0])-1 && aSort[i]!=0; i++){
      pBox = mergeComptBox(aSort[i], pBox);
      aSort[i] = 0;
    }
    aSort[i] = mergeComptBox(aSort[i], pBox);
  }
  for(i=0; i<sizeof(aSort)/sizeof(aSort[0]); i++){
    p->pComptBox = mergeComptBox(aSort[i], p->pComptBox);
  }
}

/*
** Test to see if the point x,y is contained within the given
** boundary or is on the outside of the boundary.
*/
static int pointWithinBoundary(
  Boundary *aBound,      /* The boundary */
  int nBound,            /* Number of segments in the boundary */
  double x, double y     /* The point to test */
){
  int inside = 0;
  int i;
  for(i=0; i<nBound; i++){
    double x0, y0, x1, y1;
    Segment *p = aBound[i].pSeg;
    x0 = p->from[X_IDX];
    y0 = p->from[Y_IDX];
    x1 = p->to[X_IDX];
    y1 = p->to[Y_IDX];
    if( x0==x1 ) continue;
    if( (x0>x && x1>x) || (x0<x && x1<x) ) continue;
    if( y1 - (x1-x)*(y1-y0)/(x1-x0) >= y ) inside = !inside;
  }
  return inside;
}

/*
** Find a boundary which contains xI, yI.  If the size of the boundary
** is set to 0, that means no such boundary exists.
*/
static int findBoundary(
  Wallset *p,             /* The wallset */
  double xI, double yI,   /* A point that the boundary should be near */
  int nBound,             /* Number of slots available in aBound[] */
  Boundary *aBound        /* OUT: Write results here */
){
  int n = 0;
  ComptBox *pBox;

  buildComptBoxCache(p);
  for(pBox=p->pComptBox; pBox; pBox=pBox->pNext){
    if( xI<pBox->bbox.l || xI>pBox->bbox.r || yI<pBox->bbox.b || yI>pBox->bbox.t ) continue;
    aBound[0] = pBox->prim;
    n = completeBoundary(p, 0, nBound, aBound);
    if( n>0 && pointWithinBoundary(aBound, n, xI, yI) ) break;
    n = 0;
  }
  return n;
}


/*
** Do an check of the integrity of the internal data structures.  If
** a problem is found, leave an error message in interp->result and
** return TCL_ERROR.  Return TCL_OK if everything is OK.
*/
static int selfCheck(Tcl_Interp *interp, Wallset *p){
  Link *pLink;
  Segment *pSeg;
  int h;
  char zErr[200];

  for(pLink=p->pAll; pLink; pLink=pLink->pNext){
    pSeg = pLink->pLinkNode;
    h = hashInt(pSeg->id);
    if(!segmentOnList(pSeg, p->hashId[h]) ){
      sprintf(zErr, "segment %d missing from hashId[%d]", pSeg->id, h);
      Tcl_SetResult(interp, zErr, TCL_VOLATILE);
      return TCL_ERROR;
    }
    h = hashCoord(pSeg->from[X_IDX], pSeg->from[Y_IDX]);
    if(!segmentOnList(pSeg, p->hashFrom[h]) ){
      sprintf(zErr, "segment %d missing from hashFrom[%d]", pSeg->id, h);
      Tcl_SetResult(interp, zErr, TCL_VOLATILE);
      return TCL_ERROR;
    }
    h = hashCoord(pSeg->to[X_IDX], pSeg->to[Y_IDX]);
    if(!segmentOnList(pSeg, p->hashTo[h]) ){
      sprintf(zErr, "segment %d missing from hashTo[%d]", pSeg->id, h);
      Tcl_SetResult(interp, zErr, TCL_VOLATILE);
      return TCL_ERROR;
    }
  }
  return TCL_OK;
}

/*
** The maximum number of segments in a boundary
*/
#define MX_BOUND 1000


/*
** This routine runs when a method is executed against a wallset.
*/
static int wallsetMethodProc(
  void *pArg,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  Wallset *p = (Wallset*)pArg;
  Boundary aBound[MX_BOUND];

#if 0
  /* For debugging....
  ** Print each wallset command before it is executed.
  */
  { int i;
    for(i=0; i<objc; i++){
      printf("%s%c", Tcl_GetStringFromObj(objv[i], 0), i<objc-1 ? ' ' : '\n');
    }
  }
#endif


  /* The following bit of magic implements a switch() statement that
  ** invokes the correct code based on the value of objv[1].  The
  ** mktclopts.tcl script scans this source file looking for "case"
  ** statements then generates the "wallset.h" include file that contains
  ** all of the code necessary to implement the switch.
  **
  ** In this way, we can add new methods to the command simply by adding
  ** new cases within the switch.  All the related switch code is
  ** regenerated automatically.
  */
  {
#include "wallset_cases.h"
  {

  /***************************************************************************
  ** Command methods
  */

  /*
  ** tclmethod:  WALLSET atvertex X Y
  ** title:   Return a list of wall IDs that connect to vertex X,Y
  */
  case WALLSET_ATVERTEX: {
    Link *pList;
    double x, y;
    Tcl_Obj *pResult;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y) ) return TCL_ERROR;
    pResult = Tcl_NewObj();
    ignoreNone(p);
    pList = segmentsAtVertex(p, x*p->rXZoom, y*p->rYZoom);
    while( pList ){
      Segment *pSeg=pList->pLinkNode;
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
      pList = pList->pNext;
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET boundary X Y
  ** title:   Return indices of segments forming a boundary around X Y
  */
  case WALLSET_BOUNDARY: {
    int nBound;
    double x, y;
    Tcl_Obj *pResult;
    int i;
    int showDetail = 0;

    if( objc==5 && strcmp(Tcl_GetStringFromObj(objv[2],0),"-detail")==0 ){
      showDetail = 1;
      objc--;
      objv++;
    }
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?-detail? X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y) ) return TCL_ERROR; 
    nBound = findBoundary(p, x*p->rXZoom, y*p->rYZoom, MX_BOUND, aBound);
    if( nBound>MX_BOUND ) nBound = 0;
    pResult = Tcl_NewObj();
    for(i=0; i<nBound; i++){
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(aBound[i].pSeg->id));
      if( showDetail ){
        Tcl_ListObjAppendElement(0, pResult,
           ODIE_CONSTANT_STRING(aBound[i].backwards ? "right" : "left"));
      }
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET closure WALLID BACKWARDS CHECKPRIMARY ?-coords?
  ** title:   Return the closure of a wall
  **
  ** The closure is a path of walls going clockwise from the wall given.
  ** The return value is a list consisting of wall IDs alternating with
  ** keywords "left" or "right" indicating which side of the wall applies.
  ** If the CHECKPRIMARY flag is true and the WALLID/BACKWARDS is not the
  ** primary wall id for the closure, then return an empty string.  The
  ** primary wall id is the wall id with the lowest id number, or if
  ** two walls in the closure have the same id, then the one that goes
  ** on the right side of the wall.
  */
  case WALLSET_CLOSURE: {
    int id;
    int nBound, i, checkPrim;
    Tcl_Obj *pResult;
    int coordsFlag = 0;
    int noerrFlag = 0;
    if( objc!=5 && objc!=6 ){
      Tcl_WrongNumArgs(interp, 2, objv, 
           "WALLID BACKWARDS CHECKPRIMARY ?-coords? ?-noerr?");
      return TCL_ERROR;
    }
    if( objc==6 ){
      const char *zOpt = Tcl_GetStringFromObj(objv[5],0);
      if( strcmp(zOpt,"-coords")==0 ){
        coordsFlag = 1;
      }else if( strcmp(zOpt,"-noerr")==0 ){
        noerrFlag = 1;
      }else{
        Tcl_AppendResult(interp, "unknown option: ", zOpt, 0);
        return TCL_ERROR;
      }
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &id) ) return TCL_ERROR;
    if( (aBound[0].pSeg = findSegment(p, id))==0 ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[2],0), " does not exist", 0);
      return TCL_ERROR;
    }
    if( Tcl_GetBooleanFromObj(interp, objv[3], &aBound[0].backwards) ){
      return TCL_ERROR;
    }
    if( Tcl_GetBooleanFromObj(interp, objv[4], &checkPrim) ){
      return TCL_ERROR;
    }
    ignoreNone(p);
    nBound = completeBoundary(p, checkPrim, MX_BOUND, aBound);  
    pResult = Tcl_NewObj();
    if( nBound<0 && noerrFlag ) nBound = -nBound;
    for(i=0; i<nBound; i++){
      if( coordsFlag ){
        double x, y;
        Segment *pSeg = aBound[i].pSeg;
        if( aBound[i].backwards ){
          x = pSeg->to[X_IDX];
          y = pSeg->to[Y_IDX];
        }else{
          x = pSeg->from[X_IDX];
          y = pSeg->from[Y_IDX];
        }
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(x/p->rXZoom));
        Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(y/p->rYZoom));
      }else{
        Tcl_ListObjAppendElement(0, pResult,
                                       Tcl_NewIntObj(aBound[i].pSeg->id));
        Tcl_ListObjAppendElement(0, pResult,
             ODIE_CONSTANT_STRING(aBound[i].backwards ? "right" : "left"));
      }
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET comptlist
  ** title:   Return a list of all compartments
  **
  ** A compartment is a closed circuit of walls.  This routine returns
  ** a list of all compartments.  Each element of the list consists of
  ** the primary wall for the compartment followed by a bounding box
  ** for the compartment.
  */
  case WALLSET_COMPTLIST: {
    ComptBox *pBox;
    Tcl_Obj *pResult = Tcl_NewObj();
    buildComptBoxCache(p);
    for(pBox=p->pComptBox; pBox; pBox=pBox->pNext){
      Tcl_Obj *pElem = Tcl_NewObj();
      Tcl_ListObjAppendElement(0, pElem,Tcl_NewIntObj(pBox->prim.pSeg->id));
      Tcl_ListObjAppendElement(0, pElem, Tcl_NewIntObj(pBox->prim.backwards));
      Tcl_ListObjAppendElement(0, pElem, Tcl_NewDoubleObj(pBox->bbox.l/p->rXZoom));
      Tcl_ListObjAppendElement(0, pElem, Tcl_NewDoubleObj(pBox->bbox.b/p->rYZoom));
      Tcl_ListObjAppendElement(0, pElem, Tcl_NewDoubleObj(pBox->bbox.r/p->rXZoom));
      Tcl_ListObjAppendElement(0, pElem, Tcl_NewDoubleObj(pBox->bbox.t/p->rYZoom));
      Tcl_ListObjAppendElement(0, pResult, pElem);
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET primary X Y
  ** title:   Return the primary segment of the compartment enclosing X,Y
  **
  ** The primary segment is the segment with the smallest ID.  If the
  ** same segment occurs twice on the list (in other words, if the
  ** same compartment is on both sides of a wall), then the right side
  ** (as measured facing the direction of travel from x0,y0 -> x1,y1) 
  ** is used.
  */
  case WALLSET_PRIMARY: {
    int nBound;
    double x, y;
    int i, sideSmallest;
    int idSmallest;
    Tcl_Obj *pResult;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y) ) return TCL_ERROR; 
    nBound = findBoundary(p, x*p->rXZoom, y*p->rYZoom, MX_BOUND, aBound);
    if( nBound>0 && nBound<MX_BOUND ){
      idSmallest = aBound[0].pSeg->id;
      sideSmallest = aBound[0].backwards;
      for(i=1; i<nBound; i++){
        if( aBound[i].pSeg->id>idSmallest ) continue;
        if( aBound[i].pSeg->id<idSmallest || !sideSmallest ){
          idSmallest = aBound[i].pSeg->id;
          sideSmallest = aBound[i].backwards;
        }
      }
      pResult = Tcl_NewObj();
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(idSmallest));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(sideSmallest));
      Tcl_SetObjResult(interp, pResult);
    }
    break;
  }

  /*
  ** tclmethod:  WALLSET corners X Y
  ** title:   Return vertices of compartment containing X,Y
  */
  case WALLSET_CORNERS: {
    int nBound, i;
    double x, y;
    Tcl_Obj *pResult;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y) ) return TCL_ERROR; 
    nBound = findBoundary(p, x*p->rXZoom, y*p->rYZoom, MX_BOUND, aBound);
    if( nBound>MX_BOUND ) nBound = 0;
    pResult = Tcl_NewObj();
    for(i=0; i<nBound; i++){
      Segment *pSeg = aBound[i].pSeg;
      if( aBound[i].backwards ){
        x = pSeg->to[X_IDX];
        y = pSeg->to[Y_IDX];
      }else{
        x = pSeg->from[X_IDX];
        y = pSeg->from[Y_IDX];
      }
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(x/p->rXZoom));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(y/p->rYZoom));
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET delete ID
  ** title:   Delete a single segment of a wall given by ID
  */
  case WALLSET_DELETE: {
    int id;
    Segment *pSeg;
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "ID");
      return TCL_ERROR;
    }
    if( p->busy ){
      Tcl_AppendResult(interp, "cannot \"delete\" from within a \"foreach\"",0);
      return TCL_ERROR;
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &id) ) return TCL_ERROR;
    if( (pSeg = findSegment(p, id))==0 ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[2],0), " does not exist", 0);
      return TCL_ERROR;
    }
    clearComptBoxCache(p);
    LinkRemove(&pSeg->pAll);
    /* We intentionally do not remove pSeg->pSet because it might not be
    ** a well-formed list */
    LinkRemove(&pSeg->pHash);
    LinkRemove(&pSeg->pFrom);
    LinkRemove(&pSeg->pTo);
    Odie_Free((char *)pSeg);
    break;
  }

  /*
  ** tclmethod:  WALLSET destroy
  ** title:   Destroy this wallset
  */
  case WALLSET_DESTROY: {
    Tcl_DeleteCommand(interp,Tcl_GetString(objv[0]));
    break;
  }
  
  /*
  ** tclmethod:  WALLSET firstboundary X Y
  ** title:   Find a wall on the boundary of compartment containing X Y
  **
  ** Returns a list of two elements on success or an empty list if no
  ** suitable boundary could be found.  The first element is the ID of a
  ** wall that forms part of the boundary for the compartment containing
  ** point X,Y.  The second element is TRUE if X,Y is to the right of the
  ** wall and false if it is to the left.
  **
  ** The right/left designation assumes a right-handed coordinate system.
  */
  case WALLSET_FIRSTBOUNDARY: {
    int isBack;
    Segment *pSeg;
    double x, y;
    int rc;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y) ) return TCL_ERROR;
    ignoreNone(p);
    rc = firstBoundarySegment(p, x*p->rXZoom, y*p->rYZoom, &pSeg, &isBack);
    if( rc==0 ){
      Tcl_Obj *pResult = Tcl_NewObj();
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(isBack));
      Tcl_SetObjResult(interp, pResult);
    }
    break;
  }

  /*
  ** tclmethod: WALLSET foreach CODE
  ** title:  Run CODE for each segment of the wallset
  */
  case WALLSET_FOREACH: {
    Link *pLink;
    int rc = TCL_OK;
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "CODE");
      return TCL_ERROR;
    }
    p->busy++;
    for(pLink=p->pAll; pLink && rc==TCL_OK; pLink=pLink->pNext){
      Segment *pSeg = pLink->pLinkNode;
      Tcl_SetVar2Ex(interp, "x0", 0, Tcl_NewDoubleObj(pSeg->from[X_IDX]/p->rXZoom), 0);
      Tcl_SetVar2Ex(interp, "y0", 0, Tcl_NewDoubleObj(pSeg->from[Y_IDX]/p->rYZoom), 0);
      Tcl_SetVar2Ex(interp, "x1", 0, Tcl_NewDoubleObj(pSeg->to[X_IDX]/p->rXZoom), 0);
      Tcl_SetVar2Ex(interp, "y1", 0, Tcl_NewDoubleObj(pSeg->to[Y_IDX]/p->rYZoom), 0);
      Tcl_SetVar2Ex(interp, "id", 0, Tcl_NewIntObj(pSeg->id), 0);
      Tcl_SetVar2Ex(interp, "lc", 0, Tcl_NewIntObj(pSeg->idLC), 0);
      Tcl_SetVar2Ex(interp, "rc", 0, Tcl_NewIntObj(pSeg->idRC), 0);
      Tcl_SetVar2Ex(interp, "virtual", 0, Tcl_NewIntObj(pSeg->isBoundary), 0);
      rc = Tcl_EvalObjEx(interp, objv[2], 0);
    }
    if( rc==TCL_BREAK ) rc = TCL_OK;
    p->busy--;
    return rc;
  }

  /*
  ** tclmethod: WALLSET info ID
  ** title:  Return information about a single wall segment
  */
  case WALLSET_INFO: {
    int id;
    Segment *pSeg;
    Tcl_Obj *pResult;
    if( objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "ID");
      return TCL_ERROR;
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &id) ) return TCL_ERROR;
    if( (pSeg = findSegment(p, id))==0 ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[2],0), " does not exist", 0);
      return TCL_ERROR;
    }
    pResult = Tcl_NewObj();
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(pSeg->from[X_IDX]/p->rXZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(pSeg->from[Y_IDX]/p->rYZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(pSeg->to[X_IDX]/p->rXZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(pSeg->to[Y_IDX]/p->rYZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->idLC));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->idRC));
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /* tclmethod:  WALLSET insert X0 Y0 X1 Y1 ID LC RC VIRTUAL
  ** title:   Create a new wall within the wallset
  */
  case WALLSET_INSERT: {
    int id;
    int h,virtual=0;
    double x0, y0, x1, y1;
    int idLC, idRC;
    Segment *pSeg;
    if( objc!=9 && objc!=10){
      Tcl_WrongNumArgs(interp, 2, objv, "X0 Y0 X1 Y1 ID LC RC ?1|0?");
      return TCL_ERROR;
    }
    if( p->busy ){
      Tcl_AppendResult(interp, "cannot \"insert\" from within a \"foreach\"",0);
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &x1) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[5], &y1) ) return TCL_ERROR;
    if( Tcl_GetIntFromObj(interp, objv[6], &id) ) return TCL_ERROR;
    if( Tcl_GetIntFromObj(interp, objv[7], &idLC) ) return TCL_ERROR;
    if( Tcl_GetIntFromObj(interp, objv[8], &idRC) ) return TCL_ERROR;
    if(objc==10) {
      if( Tcl_GetIntFromObj(interp, objv[8], &virtual) ) return TCL_ERROR;
    }
    x0 = roundCoord(x0*p->rXZoom);
    y0 = roundCoord(y0*p->rYZoom);
    x1 = roundCoord(x1*p->rXZoom);
    y1 = roundCoord(y1*p->rYZoom);
    if( findSegment(p, id) ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[6],0), " already exists", 0);
      return TCL_ERROR;
    }
    if( floatCompare(x0,x1)==0 && floatCompare(y0,y1)==0 ){
      /* Tcl_AppendResult(interp, "endpoints must be distinct", 0); */
      /* return TCL_ERROR; */
      return TCL_OK;  /* Not an error.  Just a no-op. */
    }
    clearComptBoxCache(p);
    pSeg = (Segment *)Odie_Alloc( sizeof(*pSeg) );
    if( pSeg==0 ) return TCL_ERROR;
    pSeg->id = id;
    pSeg->idLC = idLC;
    pSeg->idRC = idRC;
    pSeg->from[X_IDX] = x0;
    pSeg->from[Y_IDX] = y0;
    pSeg->to[X_IDX] = x1;
    pSeg->to[Y_IDX] = y1;
    pSeg->isBoundary=virtual;

    LinkInit(pSeg->pAll, pSeg);
    LinkInit(pSeg->pSet, pSeg);
    LinkInit(pSeg->pHash, pSeg);
    LinkInit(pSeg->pFrom, pSeg);
    LinkInit(pSeg->pTo, pSeg);
    LinkInsert(&p->pAll, &pSeg->pAll);
    h = hashInt(id);
    LinkInsert(&p->hashId[h], &pSeg->pHash);
    h = hashCoord(pSeg->from[X_IDX], pSeg->from[Y_IDX]);
    LinkInsert(&p->hashFrom[h], &pSeg->pFrom);
    h = hashCoord(pSeg->to[X_IDX], pSeg->to[Y_IDX]);
    LinkInsert(&p->hashTo[h], &pSeg->pTo);
    break;
  }

  /*
  ** tclmethod:  WALLSET intersect X0 Y0 X1 Y1
  ** title:   Find the intersection of X0,Y0->X1,Y1 with a segment
  **
  ** Scan all segments in the wallset looking for one that intersects with
  ** a line from X0,Y0 to X1,Y1.  If the intersection occurs at x0,y0, it
  ** is ignored, but intersections at x1,y1 count.  If no such intersection
  ** exists, return the empty string.  If there are one or more intersections,
  ** return the ID of the segment and the X and Y coordinates of the nearest
  ** intersection to X0,Y0.  
  */
  case WALLSET_INTERSECT: {
    double x0,y0,x1,y1;
    double adx, ady;
    Link *pI;
    int id;
    double nx, ny;
    double mindist2 = -1.0;
    if( objc!=6 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X0 Y0 X1 Y1");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &x1) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[5], &y1) ) return TCL_ERROR;
    x0 = roundCoord(x0*p->rXZoom);
    y0 = roundCoord(y0*p->rYZoom);
    x1 = roundCoord(x1*p->rXZoom);
    y1 = roundCoord(y1*p->rYZoom);
    adx = x1-x0;
    ady = y1-y0;
    if( adx==0.0 && ady==0.0 ) break;
    for(pI=p->pAll; pI; pI=pI->pNext){
      double bdx, bdy, denom, num1;
      Segment *pSeg;
      pSeg = pI->pLinkNode;
      bdx = pSeg->to[X_IDX] - pSeg->from[X_IDX];
      bdy = pSeg->to[Y_IDX] - pSeg->from[Y_IDX];
      denom = adx*bdy - ady*bdx;
      num1 = (y0-pSeg->from[Y_IDX])*bdx - (x0-pSeg->from[X_IDX])*bdy;
      if( denom==0.0 ){
        /* The reference line and segment are parallel */
        if( num1==0.0 ){
          /* The reference line and segment are colinear */
          if( samePoint(x0,y0,pSeg->from[X_IDX],pSeg->from[Y_IDX])
                  && adx*bdx<=0.0 && ady*bdy<=0.0 ){
            continue;
          }
          if( samePoint(x0,y0,pSeg->to[X_IDX],pSeg->to[Y_IDX])
                  && adx*bdx>=0.0 && ady*bdy>=0.0 ){
            continue;
          }
          if( between(pSeg->from[Y_IDX],y0,y1) && between(pSeg->from[X_IDX],x0,x1) ){
            double dx, dy, dist2;
            dx = pSeg->from[X_IDX] - x0;
            dy = pSeg->from[Y_IDX] - y0;
            dist2 = dx*dx + dy*dy;
            if( mindist2<0 || mindist2>dist2 ){
              mindist2 = dist2;
              nx = pSeg->from[X_IDX];
              ny = pSeg->from[Y_IDX];
              id = pSeg->id;
            }
          }
          if( between(pSeg->to[Y_IDX],y0,y1) && between(pSeg->to[X_IDX],x0,x1) ){
            double dx, dy, dist2;
            dx = pSeg->to[X_IDX] - x0;
            dy = pSeg->to[Y_IDX] - y0;
            dist2 = dx*dx + dy*dy;
            if( mindist2<0 || mindist2>dist2 ){
              mindist2 = dist2;
              nx = pSeg->to[X_IDX];
              ny = pSeg->to[Y_IDX];
              id = pSeg->id;
            }
          }
          if( between(y0,pSeg->from[Y_IDX],pSeg->to[Y_IDX]) && between(x0,pSeg->from[X_IDX],pSeg->to[X_IDX]) ){
            if( mindist2<0 || mindist2>0.0 ){
              mindist2 = 0.0;
              nx = x0;
              ny = y0;
              id = pSeg->id;
            }
          }
        }
      }else{
        /* The reference line and segment are not parallel */
        double r, s;
        r = num1/denom;
        s = ((y0-pSeg->from[Y_IDX])*adx - (x0-pSeg->from[X_IDX])*ady)/denom;
        if( r>0 && r<=1.0 && s>=0.0 && s<=1.0 ){
          double dx, dy, dist2;
          dx = r*adx;
          dy = r*ady;
          dist2 = dx*dx + dy*dy;
          if( dist2>=GRAIN && (mindist2<0 || mindist2>dist2) ){
            mindist2 = dist2;
            nx = x0 + dx;
            ny = y0 + dy;
            id = pSeg->id;
          }
        }
      }
    }
    if( mindist2>=0.0 ){
      Tcl_Obj *pResult;
      pResult = Tcl_NewObj();
      nx = roundCoord(nx);
      ny = roundCoord(ny);
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(nx/p->rXZoom));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(ny/p->rYZoom));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(id));
      Tcl_SetObjResult(interp, pResult);
    }
    break;
  }

  /*
  ** tclmethod:  WALLSET left ID LC
  ** title:   Change the left compartment of a line segment
  */
  case WALLSET_LEFT: {
    int id, idLC;
    Segment *pSeg;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "ID LC");
      return TCL_ERROR;
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &id) ) return TCL_ERROR;
    if( Tcl_GetIntFromObj(interp, objv[3], &idLC) ) return TCL_ERROR;
    if( (pSeg = findSegment(p, id))==0 ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[2],0), " does not exist", 0);
      return TCL_ERROR;
    }
    pSeg->idLC = idLC;
    break;
  }

  /*
  ** tclmethod: WALLSET list
  ** title:  Return a list of all wall segment identifiers
  */
  case WALLSET_LIST: {
    Link *pLink;
    Tcl_Obj *pResult;
    pResult = Tcl_NewObj();
    for(pLink=p->pAll; pLink; pLink=pLink->pNext){
      Segment *pSeg=pLink->pLinkNode;
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET looseends
  ** title:   Return a list of walls that are have unconnected ends
  **
  ** For each unconnected end, the list contains four elements:
  **    1.  The wallid
  **    2.  0 for the "from" end, "1" for the "to" end
  **    3.  The X coordinate of the loose end
  **    4.  The Y coordinate of the loose end
  */
  case WALLSET_LOOSEENDS: {
    Segment *pSeg;
    Link *pAll, *pList;
    Tcl_Obj *pRes = Tcl_NewObj();
    for(pAll=p->pAll; pAll; pAll=pAll->pNext){
      pSeg = pAll->pLinkNode;
      pList = segmentsAtVertex(p, pSeg->from[X_IDX], pSeg->from[Y_IDX]);
      if( LinkCount(pList)==1 ){
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewIntObj(pSeg->id));
        Tcl_ListObjAppendElement(0, pRes, ODIE_INT_ZERO());
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewDoubleObj(pSeg->from[X_IDX]/p->rXZoom));
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewDoubleObj(pSeg->from[Y_IDX]/p->rYZoom));
      }
      pList = segmentsAtVertex(p, pSeg->to[X_IDX], pSeg->to[Y_IDX]);
      if( LinkCount(pList)==1 ){
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewIntObj(pSeg->id));
        Tcl_ListObjAppendElement(0, pRes, ODIE_INT_ONE());
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewDoubleObj(pSeg->to[X_IDX]/p->rXZoom));
        Tcl_ListObjAppendElement(0, pRes, Tcl_NewDoubleObj(pSeg->to[Y_IDX]/p->rYZoom));
      }
    }
    Tcl_SetObjResult(interp, pRes);
    break;
  }

  /*
  ** tclmethod:  WALLSET nearest vertex|point X Y
  ** title:   Find the nearest vertex or point to a point in the plan
  */
  case WALLSET_NEAREST: {
    int type;
    double x, y, near_x, near_y;
    static const char *NEAR_strs[] = { "point", "vertex", 0 };
    enum NEAR_enum { NEAR_POINT, NEAR_VERTEX, };
    Link *pLink;
    Tcl_Obj *pResult;
    double dx, dy, dist;

    if( objc!=5 ){
      Tcl_WrongNumArgs(interp, 2, objv, "point|vertex X Y");
      return TCL_ERROR;
    }
    if( Tcl_GetIndexFromObj(interp, objv[2], NEAR_strs, "option", 0, &type) ){
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[3], &x) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &y) ) return TCL_ERROR;
    x *= p->rXZoom;
    y *= p->rYZoom;
    ignoreNone(p);
    if( type==NEAR_POINT ){
      pLink = nearestPoint(p, x, y, &near_x, &near_y);
    }else if( type==NEAR_VERTEX ){
      pLink = nearestVertex(p, x, y, &near_x, &near_y);
    }else{
      /* Cannot happen */ return TCL_ERROR;
    }
    if( pLink==0 ) return TCL_OK;  /* There are not segments in the wallset */
    pResult = Tcl_NewObj();
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(near_x/p->rXZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(near_y/p->rYZoom));
    dx = x - near_x;
    dy = y - near_y;
    dist = sqrt(dx*dx + dy*dy);
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewDoubleObj(dist/p->rXZoom));
    Tcl_ListObjAppendElement(0, pResult, Tcl_NewObj());
    while( pLink ){
      Segment *pSeg=pLink->pLinkNode;
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
      pLink = pLink->pNext;
    }
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /*
  ** tclmethod:  WALLSET nextcwwall X0 Y0 X1 Y1
  ** title:   Find a wall on X1,Y1 clockwise from X0,Y0->X1,Y1
  */
  case WALLSET_NEXTCWWALL: {
    int isBack;
    Segment *pSeg;
    double x0, y0, x1, y1;
    int rc;
    if( objc!=6 ){
      Tcl_WrongNumArgs(interp, 2, objv, "X0 Y0 X1 Y1");
      return TCL_ERROR;
    }
    if( Tcl_GetDoubleFromObj(interp, objv[2], &x0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[3], &y0) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[4], &x1) ) return TCL_ERROR;
    if( Tcl_GetDoubleFromObj(interp, objv[5], &y1) ) return TCL_ERROR;
    x0 = roundCoord(x0*p->rXZoom);
    y0 = roundCoord(y0*p->rYZoom);
    x1 = roundCoord(x1*p->rXZoom);
    y1 = roundCoord(y1*p->rYZoom);
    rc = nextCwSegment(p, x0, y0, x1, y1, &pSeg, &isBack);
    if( rc==0 ){
      Tcl_Obj *pResult = Tcl_NewObj();
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(pSeg->id));
      Tcl_ListObjAppendElement(0, pResult, Tcl_NewIntObj(isBack));
      Tcl_SetObjResult(interp, pResult);
    }
    break;
  }

  /*
  ** tclmethod:  WALLSET right ID RC
  ** title:   Change the right compartment of a line segment
  */
  case WALLSET_RIGHT: {
    int id, idRC;
    Segment *pSeg;
    if( objc!=4 ){
      Tcl_WrongNumArgs(interp, 2, objv, "ID RC");
      return TCL_ERROR;
    }
    if( Tcl_GetIntFromObj(interp, objv[2], &id) ) return TCL_ERROR;
    if( Tcl_GetIntFromObj(interp, objv[3], &idRC) ) return TCL_ERROR;
    if( (pSeg = findSegment(p, id))==0 ){
      Tcl_AppendResult(interp, "segment ",
        Tcl_GetStringFromObj(objv[2],0), " does not exist", 0);
      return TCL_ERROR;
    }
    pSeg->idRC = idRC;
    break;
  }

  /*
  ** tclmethod: WALLSET segments
  ** title:  Write all wall segments out into the segset native datatype
  */
  case WALLSET_SEGMENTS: {
    
  }

  /*
  ** tclmethod:  WALLSET selfcheck
  ** title:   Verify the integrity of internal data structures
  */
  case WALLSET_SELFCHECK: {
    return selfCheck(interp, p);
  }

  /*
  ** tclmethod:  WALLSET zoom ?ZOOM?
  ** title:   Query or change the zoom factor.
  */
  case WALLSET_ZOOM: {
    Tcl_Obj *pResult;
    if( objc!=2 && objc!=3 ){
      Tcl_WrongNumArgs(interp, 2, objv, "?ZOOM?");
      return TCL_ERROR;
    }
    if( objc==3 ){
      double r;
      if( Tcl_GetDoubleFromObj(interp, objv[2], &r) ) return TCL_ERROR;
      if( r==0.0 ){
        Tcl_AppendResult(interp, "zoom must be non-zero", 0);
        return TCL_ERROR;
      }
      p->rYZoom = r;
      p->rXZoom = fabs(r);
    }
    pResult = Tcl_NewDoubleObj(p->rYZoom);
    Tcl_SetObjResult(interp, pResult);
    break;
  }

  /* End of the command methods.  The brackets that follow terminate the
  ** automatically generated switch.
  ****************************************************************************/
  }
  }

#if 0
  /* Sanity checking for debugging */
  if( selfCheck(interp, p) ){
    return TCL_ERROR;
  }
#endif
  return TCL_OK;
}

/*
** tclcmd: wallset WALLSET
** title: Create a new wallset object
** This routine runs when the "wallset" command is invoked to create a
** new wallset.
*/
int Odie_WallsetCreateProc(
  void *NotUsed,
  Tcl_Interp *interp,
  int objc,
  Tcl_Obj *CONST objv[]
){
  char *zCmd;
  Wallset *p;
  if( objc!=2 ){
    Tcl_WrongNumArgs(interp, 1, objv, "WALLSET");
    return TCL_ERROR;
  }
  zCmd = Tcl_GetStringFromObj(objv[1], 0);
  p = (Wallset *)Odie_Alloc( sizeof(*p) );
  p->rXZoom = 100.0;
  p->rYZoom = -100.0;
  Tcl_CreateObjCommand(interp, zCmd, wallsetMethodProc, p, destroyWallset);
  return TCL_OK;
}