r327: Added textfile output
[ctsim.git] / libctsim / phantom.cpp
index ec5d0dcf6e0525cc38af1f1365091dbb81889e4b..d8fa862f625724968e328fc86be4eabdb03deaea 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phantom.cpp,v 1.3 2000/06/19 20:08:09 kevin Exp $
+**  $Id: phantom.cpp,v 1.24 2001/01/02 05:33:37 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
 
 #include "ct.h"
 
+const int PhantomElement::POINTS_PER_CIRCLE = 360;\r
+const double PhantomElement::SCALE_PELEM_EXTENT=0.005;  // increase pelem limits by 0.5% \r
+
+const int Phantom::PHM_INVALID = -1;
+const int Phantom::PHM_HERMAN = 0;
+const int Phantom::PHM_B_HERMAN = 1;
+const int Phantom::PHM_SHEPP_LOGAN = 2;
+const int Phantom::PHM_B_SHEPP_LOGAN = 3;
+const int Phantom::PHM_UNITPULSE = 4;
+
+const char* Phantom::s_aszPhantomName[] = 
+{
+  {"herman"},
+  {"herman-b"},
+  {"shepp-logan"},
+  {"shepp-logan-b"},
+  {"unitpulse"},
+};
+
+const char* Phantom::s_aszPhantomTitle[] = 
+{
+  {"Herman Head"},
+  {"Herman Head (Bordered)"},
+  {"Shepp-Logan"},
+  {"Shepp-Logan (Bordered)"},
+  {"Unit Pulse"},
+};
+
+const int Phantom::s_iPhantomCount = sizeof(s_aszPhantomName) / sizeof(const char*);
+
 
 // CLASS IDENTIFICATION
-//   Phanton
+//   Phantom
 //
 
-Phantom::Phantom (void)
+Phantom::Phantom ()
+{
+  init ();
+}
+
+
+Phantom::Phantom (const char* const phmName)
+{
+  init ();
+  createFromPhantom (phmName);
+}
+
+void 
+Phantom::init ()
 {
   m_nPElem = 0;
   m_xmin = 1E30;
@@ -41,10 +84,11 @@ Phantom::Phantom (void)
   m_ymax = -1E30;
   m_diameter = 0;
   m_composition = P_PELEMS;
+  m_fail = false;
+  m_id = PHM_INVALID;
 }
 
-
-Phantom::~Phantom (void)
+Phantom::~Phantom ()
 {
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
     delete *i;
@@ -52,29 +96,91 @@ Phantom::~Phantom (void)
 }
 
 
-void
-Phantom::create (const int phmid)
+const char*
+Phantom::convertPhantomIDToName (int phmID)
+{
+  static const char *name = "";
+
+  if (phmID >= 0 && phmID < s_iPhantomCount)
+      return (s_aszPhantomName[phmID]);
+
+  return (name);
+}
+
+const char*
+Phantom::convertPhantomIDToTitle (int phmID)
+{
+  static const char *title = "";
+
+  if (phmID >= 0 && phmID < s_iPhantomCount)
+      return (s_aszPhantomName[phmID]);
+
+  return (title);
+}
+      
+int
+Phantom::convertNameToPhantomID (const char* const phmName) 
+{
+  int id = PHM_INVALID;
+
+  for (int i = 0; i < s_iPhantomCount; i++)
+      if (strcasecmp (phmName, s_aszPhantomName[i]) == 0) {
+         id = i;
+         break;
+      }
+
+  return (id);
+}
+  
+
+bool
+Phantom::createFromPhantom (const char* const phmName)
+{
+  int phmid = convertNameToPhantomID (phmName);
+  if (phmid == PHM_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid phantom name ";
+    m_failMessage += phmName;
+    return false;
+  }
+
+  m_name = phmName;
+  createFromPhantom (phmid);
+  return true;
+}
+
+bool
+Phantom::createFromPhantom (const int phmid)
 {
   switch (phmid) 
     {
-    case O_PHM_HERMAN:
+    case PHM_HERMAN:
       addStdHerman();
       break;
-    case O_PHM_ROWLAND:
-      addStdRowland();
+    case PHM_B_HERMAN:
+      addStdHermanBordered();
+      break;
+    case PHM_SHEPP_LOGAN:
+      addStdSheppLogan();
       break;
-    case O_PHM_BROWLAND:
-      addStdRowlandBordered ();
+    case PHM_B_SHEPP_LOGAN:
+      addStdSheppLoganBordered();
       break;
-    case O_PHM_UNITPULSE:
+    case PHM_UNITPULSE:
       m_composition = P_UNIT_PULSE;
       addPElem ("rectangle", 0., 0., 100., 100., 0., 0.);     // outline 
       addPElem ("ellipse", 0., 0., 1., 1., 0., 1.);          // pulse 
       break;
     default:
-      sys_error (ERR_WARNING, "Illegal phantom id %d\n", phmid);
-      break;
+      m_fail = true;
+      m_failMessage = "Illegal phantom id ";
+      m_failMessage += phmid;
+      return false;
     }
+
+  m_id = phmid;
+
+  return true;
 }
 
 
@@ -92,36 +198,43 @@ Phantom::create (const int phmid)
 bool
 Phantom::createFromFile (const char* const fname)
 {
-  bool stoploop = false;
-  bool retval = false;
+  bool bGoodFile = true;
   FILE *fp;
 
   if ((fp = fopen (fname, "r")) == NULL)
     return (false);
 
-  do {
+  m_name = fname;
+
+  while (1) {
     double cx, cy, u, v, rot, dens;
     char pelemtype[80];
-    int n = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf",
-               pelemtype, &cx, &cy, &u, &v, &rot, &dens);
+
+    int status = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf", pelemtype, &cx, &cy, &u, &v, &rot, &dens);
     
-    if (n == EOF || n == 0) {  /* end of file */
-      stoploop = true;
-      retval = false;
-    } else if (n != 7) {
-      stoploop = true;
-      retval = false;
-    } else {
-      addPElem (pelemtype, cx, cy, u, v, rot, dens);
-      retval = true;
+    if (status == static_cast<int>(EOF)) 
+      break;
+    else if (status != 7) {
+      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromFile]", fname);
+      bGoodFile = false;
     }
-  } while (stoploop == false);
+    addPElem (pelemtype, cx, cy, u, v, rot, dens);
+  }
   
   fclose (fp);
 
-  return (retval);
+  return (bGoodFile);
 }
 
+bool\r
+Phantom::fileWrite (const char* const fname)\r
+{\r
+  fstream file (fname, ios::out);\r
+\r
+  if (! file.fail())\r
+    printDefinitions (file);\r
+  return ! file.fail();\r
+}\r
 
 /* NAME
  *   addPElem          Add pelem
@@ -170,26 +283,32 @@ Phantom::addPElem (const char *type, const double cx, const double cy, const dou
  */
 
 void 
-Phantom::print (void) const
+Phantom::print (std::ostream& os) const
 {
-  printf("PRINTING Phantom\n\n");
-  printf("number of pelems in Phantom = %d\n", m_nPElem);
-  printf("limits: xmin=%8.2g  ymin=%8.2g  xmax=%8.2g  ymax=%8.2g\n",
-        m_xmin, m_ymin, m_xmax, m_ymax);
+  os << "Number of PElements: " << m_nPElem << "\n";
+  os << "Limits: xmin=" << m_xmin << ", ymin=" << m_ymin << ", xmax=" << m_xmax << ", ymax=" << m_ymax << "\n";
   
-  for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
-      printf("PELEM:\n");
-      printf("# pts=%3d        atten = %7.4f   rot = %7.2f (deg)\n",
-            (*i)->nOutlinePoints(), (*i)->atten(), convertRadiansToDegrees ((*i)->rot()));
-    
-    printf("xmin=%7.3g  ymin=%7.3g  xmax=%7.3g  ymax=%7.3g\n",
-          (*i)->xmin(), (*i)->ymin(), (*i)->xmax(), (*i)->ymax());
-    
-    //    for (int i = 0; i < m_nPoints; i++)
-    //      printf("\t%8.3g    %8.3g\n", i->xOutline()[i], i->yOutline()[i]);
+  for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {\r
+    const PhantomElement& rPE = **i;
+      os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();\r
+      os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
+      os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
+    \r
+      if (false)
+        for (int i = 0; i < rPE.nOutlinePoints(); i++)
+          os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
   }
 }
-
+\r
+void\r
+Phantom::printDefinitions (std::ostream& os) const\r
+{\r
+  for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {\r
+    const PhantomElement& rPE = **i;\r
+    rPE.printDefinition (os);\r
+  }\r
+}
+\r
 
 /* NAME
  *   show              Show vector outline of Phantom to user
@@ -200,33 +319,32 @@ Phantom::print (void) const
 
 #ifdef HAVE_SGP
 void 
-Phantom::show (void) const
+Phantom::show () const
 {
-  double wsize = m_xmax - m_xmin;
-  double xmin = m_xmin;
-  double ymin = m_ymin;
-  double xmax, ymax;
-  SGP_ID gid;
+  SGPDriver driverSGP ("Phantom Show");
+  SGP sgp (driverSGP);
 
-  if ((m_ymax - m_ymin) > wsize)
-      wsize = m_ymax - m_ymin;
-  wsize *= 1.1;
+  show (sgp);
 
-  xmax = xmin + wsize;
-  ymax = ymin + wsize; 
-  
-  printf("Drawing Phantom:\n\n");
-  printf("    data limits: %9.3g, %9.3g, %9.3g, %9.3g\n",
-        m_xmin, m_ymin, m_xmax, m_ymax);
-  printf("    window size: %9.3g, %9.3g, %9.3g, %9.3g\n",
-        xmin, ymin, xmax, ymax);
+  std::cout << "Press return to continue";
+  cio_kb_getc();
+}
 
-  gid = sgp2_init(0, 0, "Phantom Show");
-  sgp2_window (xmin, ymin, xmax, ymax);
+void 
+Phantom::show (SGP& sgp) const
+{
+  double wsize = m_xmax - m_xmin;
+  if ((m_ymax - m_ymin) > wsize) 
+      wsize = m_ymax - m_ymin;
+  wsize *= 1.01;
+  double halfWindow = wsize / 2;
 
-  draw();
+  double xcent = m_xmin + (m_xmax - m_xmin) / 2;
+  double ycent = m_ymin + (m_ymax - m_ymin) / 2;
 
-  termgrf2();
+  sgp.setWindow (xcent - halfWindow, ycent - halfWindow, xcent + halfWindow, ycent + halfWindow);
+
+  draw (sgp);
 }
 #endif
 
@@ -238,18 +356,18 @@ Phantom::show (void) const
  *   draw ()
  */
 
+#ifdef HAVE_SGP
 void 
-Phantom::draw (void) const
+Phantom::draw (SGP& sgp) const
 {
-#ifdef HAVE_SGP
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
-    sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
-#endif
+    sgp.polylineAbs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
 }
+#endif
 
 
 /* NAME
- *   addStdRowland             Make head phantom of S.W. Rowland
+ *   addStdSheppLogan  Make head phantom of Shepp-Logan
  *
  * REFERENCES
  *   S. W. Rowland, "Computer Implementation of Image Reconstruction
@@ -258,26 +376,26 @@ Phantom::draw (void) const
  */
 
 void 
-Phantom::addStdRowland (void)
+Phantom::addStdSheppLogan ()
 {
-  addPElem("ellipse",  0.0000,  0.0000, 0.6900,  0.9200,   0.0,  1.00);
-  addPElem("ellipse",  0.0000, -0.0184, 0.6624,  0.8740,   0.0, -0.98);
-  addPElem("ellipse",  0.2200,  0.0000, 0.1100,  0.3100, -18.0, -0.02);
-  addPElem("ellipse", -0.2200,  0.0000, 0.1600,  0.4100,  18.0, -0.02);
-  addPElem("ellipse",  0.0000,  0.3500, 0.2100,  0.2500,   0.0,  0.01);
-  addPElem("ellipse",  0.0000,  0.1000, 0.0460,  0.0460,   0.0,  0.01);
-  addPElem("ellipse",  0.0000, -0.1000, 0.0460,  0.0460,   0.0,  0.01);
-  addPElem("ellipse", -0.0800, -0.6050, 0.0460,  0.0230,   0.0,  0.01);
-  addPElem("ellipse",  0.0000, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
-  addPElem("ellipse",  0.0600, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
-  addPElem("ellipse",  0.5538, -0.3858, 0.0330,  0.2060, -18.0,  0.03);
+  addPElem ("ellipse",  0.0000,  0.0000, 0.6900,  0.9200,   0.0,  1.00);
+  addPElem ("ellipse",  0.0000, -0.0184, 0.6624,  0.8740,   0.0, -0.98);
+  addPElem ("ellipse",  0.2200,  0.0000, 0.1100,  0.3100, -18.0, -0.02);
+  addPElem ("ellipse", -0.2200,  0.0000, 0.1600,  0.4100,  18.0, -0.02);
+  addPElem ("ellipse",  0.0000,  0.3500, 0.2100,  0.2500,   0.0,  0.01);
+  addPElem ("ellipse",  0.0000,  0.1000, 0.0460,  0.0460,   0.0,  0.01);
+  addPElem ("ellipse",  0.0000, -0.1000, 0.0460,  0.0460,   0.0,  0.01);
+  addPElem ("ellipse", -0.0800, -0.6050, 0.0460,  0.0230,   0.0,  0.01);
+  addPElem ("ellipse",  0.0000, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
+  addPElem ("ellipse",  0.0600, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
+  addPElem ("ellipse",  0.5538, -0.3858, 0.0330,  0.2060, -18.0,  0.03);
 }
 
 void 
-Phantom::addStdRowlandBordered (void)
+Phantom::addStdSheppLoganBordered ()
 {
-  addStdRowland ();
-  addPElem ("ellipse", 0.000, 0.0000, 0.7500, 1.000, 0.0, 0.00);
+  addStdSheppLogan ();
+  addPElem ("rectangle", 0.000, 0.0000, 0.8600, 1.150, 0.0, 0.00);
 }
 
 /* NAME
@@ -289,23 +407,30 @@ Phantom::addStdRowlandBordered (void)
  */
 
 void 
-Phantom::addStdHerman (void)
-{
-  addPElem("ellipse",  0.000,  1.50,  0.375, 0.3000,  90.00, -0.003);
-  addPElem("ellipse",  0.675, -0.75,  0.225, 0.1500, 140.00,  0.010);
-  addPElem("ellipse",  0.750,  1.50,  0.375, 0.2250,  50.00,  0.003);
-  addPElem("segment",  1.375, -7.50,  1.100, 0.6250,  19.20, -0.204);
-  addPElem("segment",  1.375, -7.50,  1.100, 4.3200,  19.21,  0.204);
-  addPElem("segment",  0.000, -2.25,  1.125, 0.3750,   0.00, -0.003);
-  addPElem("segment",  0.000, -2.25,  1.125, 3.0000,   0.00,  0.003);
-  addPElem("segment", -1.000,  3.75,  1.000, 0.5000, 135.00, -0.003);
-  addPElem("segment", -1.000,  3.75,  1.000, 3.0000, 135.00,  0.003);
-  addPElem("segment",  1.000,  3.75,  1.000, 0.5000, 225.00, -0.003);
-  addPElem("segment",  1.000,  3.75,  1.000, 3.0000, 225.00,  0.003);
-  addPElem("triangle", 5.025,  3.75,  1.125, 0.5000, 110.75,  0.206);
-  addPElem("triangle",-5.025,  3.75,  1.125, 0.9000,-110.75,  0.206);
-  addPElem("ellipse",  0.000,  0.00,  8.625, 6.4687,  90.00,  0.416);
-  addPElem("ellipse",  0.000,  0.00,  7.875, 5.7187,  90.00, -0.206);
+Phantom::addStdHerman ()
+{
+  addPElem ("ellipse",  0.000,  1.50,  0.375, 0.3000,  90.00, -0.003);
+  addPElem ("ellipse",  0.675, -0.75,  0.225, 0.1500, 140.00,  0.010);
+  addPElem ("ellipse",  0.750,  1.50,  0.375, 0.2250,  50.00,  0.003);
+  addPElem ("segment",  1.375, -7.50,  1.100, 0.6250,  19.20, -0.204);
+  addPElem ("segment",  1.375, -7.50,  1.100, 4.3200,  19.21,  0.204);
+  addPElem ("segment",  0.000, -2.25,  1.125, 0.3750,   0.00, -0.003);
+  addPElem ("segment",  0.000, -2.25,  1.125, 3.0000,   0.00,  0.003);
+  addPElem ("segment", -1.000,  3.75,  1.000, 0.5000, 135.00, -0.003);
+  addPElem ("segment", -1.000,  3.75,  1.000, 3.0000, 135.00,  0.003);
+  addPElem ("segment",  1.000,  3.75,  1.000, 0.5000, 225.00, -0.003);
+  addPElem ("segment",  1.000,  3.75,  1.000, 3.0000, 225.00,  0.003);
+  addPElem ("triangle", 5.025,  3.75,  1.125, 0.5000, 110.75,  0.206);
+  addPElem ("triangle",-5.025,  3.75,  1.125, 0.9000,-110.75,  0.206);
+  addPElem ("ellipse",  0.000,  0.00,  8.625, 6.4687,  90.00,  0.416);
+  addPElem ("ellipse",  0.000,  0.00,  7.875, 5.7187,  90.00, -0.206);
+}
+
+void
+Phantom::addStdHermanBordered ()
+{
+  addStdHerman();
+  addPElem ("rectangle",  0.000, 0.00, 10.780, 8.110,  90.00, 0.000);
 }
 
 
@@ -323,11 +448,11 @@ Phantom::addStdHerman (void)
 void
 Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const
 {
-  convertToImagefile (im, in_nsample, trace, 0, im.nx());
+  convertToImagefile (im, in_nsample, trace, 0, im.nx(), true);
 }
 
 void 
-Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount) const
+Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount, bool bStoreAtColumnPos) const
 {
   int nx = im.nx();
   int ny = im.ny();
@@ -368,28 +493,45 @@ Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trac
 
   ImageFileArray v = im.getArray();
 
+  for (int ix = 0; ix < colCount; ix++)
+    for (int iy = 0; iy < ny; iy++) {
+      int iColStore = ix;
+      if (bStoreAtColumnPos)
+               iColStore += colStart;
+      v[iColStore][iy] = 0;
+    }
+
   double x_start = xmin + (colStart * xinc);
   for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
-    double x, y, xi, yi;
-    int ix, iy, kx, ky;
-    for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) {
-      for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
-       for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
-         for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
-           if ((*pelem)->isPointInside (xi, yi, PHM_COORD) == TRUE)
-             v[ix][iy] += (*pelem)->atten();
-       } // for kx
-      } /* for iy */
-    }  /* for ix */
+      const PhantomElement& rPElem = **pelem;
+      double x, y, xi, yi;
+      int ix, iy, kx, ky;
+      for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) {
+       int iColStore = ix;
+       if (bStoreAtColumnPos)
+         iColStore += colStart;
+       for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
+         for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
+           for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
+             if (rPElem.isPointInside (xi, yi, PHM_COORD) == TRUE)
+                       v[iColStore][iy] += rPElem.atten();
+             } // for kx
+         } /* for iy */
+      }  /* for ix */
   }  /* for pelem */
   
 
   if (nsample > 1) {
-    double factor = 1.0 / (nsample * nsample);
-
-    for (int ix = 0; ix < colCount; ix++)
-      for (int iy = 0; iy < ny; iy++)
-       v[ix][iy] *= factor;
+    double factor = 1.0 / static_cast<double>(nsample * nsample);
+
+\r
+    for (int ix = 0; ix < colCount; ix++) {
+               int iColStore = ix;\r
+               if (bStoreAtColumnPos)\r
+                       iColStore += colStart;\r
+               for (int iy = 0; iy < ny; iy++)
+                       v[iColStore][iy] *= factor;\r
+       }
   }
 }
 
@@ -408,20 +550,7 @@ PhantomElement::PhantomElement (const char *type, const double cx, const double
 {
   m_rot = convertDegreesToRadians (rot);   // convert angle to radians
 
-  if (strcasecmp (type, "rectangle") == 0)
-    m_type = PELEM_RECTANGLE;
-  else if (strcasecmp (type, "triangle") == 0)
-    m_type = PELEM_TRIANGLE;
-  else if (strcasecmp (type, "ellipse") == 0)
-    m_type = PELEM_ELLIPSE;
-  else if (strcasecmp (type, "sector") == 0)
-    m_type = PELEM_SECTOR;
-  else if (strcasecmp (type, "segment") == 0)
-    m_type = PELEM_SEGMENT;
-  else {
-    sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
-    m_type = PELEM_INVALID;
-  }
+  m_type = convertNameToType (type);
 
   makeTransformMatrices ();     // calc transform matrices between phantom and normalized phantomelement
   makeVectorOutline ();                // calculate vector outline of pelem 
@@ -439,15 +568,64 @@ PhantomElement::PhantomElement (const char *type, const double cx, const double
   m_rectLimits[2] = m_xmax;   m_rectLimits[3] = m_ymax;
 }
 
+\r
 
-PhantomElement::~PhantomElement (void)
+PhantomElement::~PhantomElement ()
 {
     delete m_xOutline;
     delete m_yOutline;
 }
+\r
+void\r
+PhantomElement::printDefinition (std::ostream& os) const\r
+{\r
+  os << convertTypeToName (m_type) << " " << m_cx << " " << m_cy << " " << m_u << " "\r
+    << m_v << " " << convertRadiansToDegrees (m_rot) << " " << m_atten << "\n";\r
+}\r
+
+PhmElemType
+PhantomElement::convertNameToType (const char* const typeName)
+{
+    PhmElemType type = PELEM_INVALID;
+
+    if (strcasecmp (typeName, "rectangle") == 0)
+       type = PELEM_RECTANGLE;
+    else if (strcasecmp (typeName, "triangle") == 0)
+       type = PELEM_TRIANGLE;
+    else if (strcasecmp (typeName, "ellipse") == 0)
+       type = PELEM_ELLIPSE;
+    else if (strcasecmp (typeName, "sector") == 0)
+       type = PELEM_SECTOR;
+    else if (strcasecmp (typeName, "segment") == 0)
+      type = PELEM_SEGMENT;
+    else
+       sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
+
+    return (type);
+}
+\r
+const char* const\r
+PhantomElement::convertTypeToName (PhmElemType iType)\r
+{\r
+  static char* pszType = "Unknown";\r
+\r
+  if (iType == PELEM_RECTANGLE)\r
+    pszType = "rectangle";\r
+  else if (iType == PELEM_TRIANGLE)\r
+    pszType = "triangle";\r
+  else if (iType == PELEM_ELLIPSE)\r
+    pszType = "ellipse";\r
+  else if (iType == PELEM_SECTOR)\r
+    pszType = "sector";\r
+  else if (iType == PELEM_SEGMENT)\r
+    pszType = "segment";\r
+\r
+  return pszType;\r
+}\r
+\r
 
 void 
-PhantomElement::makeTransformMatrices (void)
+PhantomElement::makeTransformMatrices ()
 {
   GRFMTX_2D temp;
 
@@ -486,10 +664,8 @@ PhantomElement::makeTransformMatrices (void)
  *   Called by phm_add_pelem()
  */
 
-static const double SCALE_PELEM_EXTENT=0.005;          // increase pelem limits by 0.5% 
-
 void
-PhantomElement::makeVectorOutline (void)
+PhantomElement::makeVectorOutline ()
 {
   double radius, theta, start, stop;
   double xfact, yfact;
@@ -578,7 +754,6 @@ PhantomElement::makeVectorOutline (void)
 }
 
 
-
 /* NAME
  *   calc_arc                  Calculate outline of a arc of a circle
  *
@@ -642,10 +817,9 @@ PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const
  */
 
 int 
-PhantomElement::numCirclePoints (double theta) const
+PhantomElement::numCirclePoints (double theta)
 {
-    if (theta < 0.0 || theta > TWOPI)
-      sys_error(ERR_WARNING, "illegal values sent to circle_pts");
+    theta = clamp (theta, 0., TWOPI);
 
     return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
 }
@@ -735,7 +909,7 @@ PhantomElement::clipLineNormalizedCoords (double& x1, double& y1, double& x2, do
 //    false if point lies outside of pelem
 
 bool
-PhantomElement::isPointInside (double x, double y, const CoordType coord_type)
+PhantomElement::isPointInside (double x, double y, const CoordType coord_type) const
 {
   if (coord_type == PHM_COORD) {
     xform_mtx2 (m_xformPhmToObj, x, y);