r589: Added threaded rasterizer
[ctsim.git] / libctsim / phantom.cpp
index a0ddc408937b8f8f382c36ae95bbc320a522ace0..dff1d0ee6bcc137ab1257d31d828921ffc937f96 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: phantom.cpp,v 1.29 2001/02/09 14:34:16 kevin Exp $
+**  $Id: phantom.cpp,v 1.30 2001/02/27 03:59:30 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
@@ -93,10 +93,10 @@ const char*
 Phantom::convertPhantomIDToName (int phmID)
 {
   static const char *name = "";
-
+  
   if (phmID >= 0 && phmID < s_iPhantomCount)
-      return (s_aszPhantomName[phmID]);
-
+    return (s_aszPhantomName[phmID]);
+  
   return (name);
 }
 
@@ -104,27 +104,27 @@ const char*
 Phantom::convertPhantomIDToTitle (int phmID)
 {
   static const char *title = "";
-
+  
   if (phmID >= 0 && phmID < s_iPhantomCount)
-      return (s_aszPhantomName[phmID]);
-
+    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);
+    if (strcasecmp (phmName, s_aszPhantomName[i]) == 0) {
+      id = i;
+      break;
+    }
+    
+    return (id);
 }
-  
+
 
 bool
 Phantom::createFromPhantom (const char* const phmName)
@@ -136,7 +136,7 @@ Phantom::createFromPhantom (const char* const phmName)
     m_failMessage += phmName;
     return false;
   }
-
+  
   m_name = phmName;
   createFromPhantom (phmid);
   return true;
@@ -146,57 +146,57 @@ bool
 Phantom::createFromPhantom (const int phmid)
 {
   switch (phmid) 
-    {
-    case PHM_HERMAN:
-      addStdHerman();
-      break;
-    case PHM_SHEPP_LOGAN:
-      addStdSheppLogan();
-      break;
-    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:
-      m_fail = true;
-      m_failMessage = "Illegal phantom id ";
-      m_failMessage += phmid;
-      return false;
-    }
-
+  {
+  case PHM_HERMAN:
+    addStdHerman();
+    break;
+  case PHM_SHEPP_LOGAN:
+    addStdSheppLogan();
+    break;
+  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:
+    m_fail = true;
+    m_failMessage = "Illegal phantom id ";
+    m_failMessage += phmid;
+    return false;
+  }
+  
   m_id = phmid;
-
+  
   return true;
 }
 
 
 /* METHOD IDENTIFICATION
- *   createFromFile          Add PhantomElements from file
- *
- * SYNOPSIS
- *   createFromFile (filename)
- *
- * RETURNS
- *   true if pelem were added
- *   false if an pelem not added
- */
+*   createFromFile          Add PhantomElements from file
+*
+* SYNOPSIS
+*   createFromFile (filename)
+*
+* RETURNS
+*   true if pelem were added
+*   false if an pelem not added
+*/
 
 bool
 Phantom::createFromFile (const char* const fname)
 {
   bool bGoodFile = true;
   FILE *fp;
-
+  
   if ((fp = fopen (fname, "r")) == NULL)
     return (false);
-
+  
   m_name = fname;
-
+  
   while (1) {
     double cx, cy, u, v, rot, dens;
     char pelemtype[80];
-
+    
     int status = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf", pelemtype, &cx, &cy, &u, &v, &rot, &dens);
     
     if (status == static_cast<int>(EOF)) 
@@ -209,7 +209,7 @@ Phantom::createFromFile (const char* const fname)
   }
   
   fclose (fp);
-
+  
   return (bGoodFile);
 }
 
@@ -217,37 +217,37 @@ bool
 Phantom::fileWrite (const char* const fname)
 {
   fstream file (fname, ios::out);
-
+  
   if (! file.fail())
     printDefinitions (file);
   return ! file.fail();
 }
 
 /* NAME
- *   addPElem          Add pelem
- *
- * SYNOPSIS
- *   addPElem (type, cx, cy, u, v, rot, atten)
- *   char *type                type of pelem (box, ellipse, etc)
- *   double cx, cy     pelem center
- *   double u,v                pelem size
- *   double rot                rotation angle of pelem (in degrees)
- *   double atten      x-ray attenuation cooefficient
- */
+*   addPElem           Add pelem
+*
+* SYNOPSIS
+*   addPElem (type, cx, cy, u, v, rot, atten)
+*   char *type         type of pelem (box, ellipse, etc)
+*   double cx, cy      pelem center
+*   double u,v         pelem size
+*   double rot         rotation angle of pelem (in degrees)
+*   double atten       x-ray attenuation cooefficient
+*/
 
 void 
 Phantom::addPElem (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
 {
   PhantomElement *pelem = new PhantomElement (type, cx, cy, u, v, rot, atten);
-
+  
   m_listPElem.push_front (pelem);
-
+  
   // update phantom limits
   if (m_xmin > pelem->xmin())    m_xmin = pelem->xmin();
   if (m_xmax < pelem->xmax())    m_xmax = pelem->xmax();
   if (m_ymin > pelem->ymin())    m_ymin = pelem->ymin();
   if (m_ymax < pelem->ymax())    m_ymax = pelem->ymax();
-
+  
   m_nPElem++;
 }
 
@@ -258,11 +258,11 @@ Phantom::addPElem (const char *type, const double cx, const double cy, const dou
 
 
 /* NAME
- *   print                             Print vertices of Phantom pelems
- *
- * SYNOPSIS
- *   print (phm)
- */
+*   print                              Print vertices of Phantom pelems
+*
+* SYNOPSIS
+*   print (phm)
+*/
 
 void 
 Phantom::print (std::ostream& os) const
@@ -272,13 +272,13 @@ Phantom::print (std::ostream& os) const
   
   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
     const PhantomElement& rPE = **i;
-      os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
-      os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
-      os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
+    os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
+    os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
+    os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
     
-      if (false)
-        for (int i = 0; i < rPE.nOutlinePoints(); i++)
-          os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
+    if (false)
+      for (int i = 0; i < rPE.nOutlinePoints(); i++)
+        os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
   }
 }
 void 
@@ -289,13 +289,13 @@ Phantom::print (std::ostringstream& os) const
   
   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
     const PhantomElement& rPE = **i;
-      os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
-      os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
-      os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
+    os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
+    os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
+    os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
     
-      if (false)
-        for (int i = 0; i < rPE.nOutlinePoints(); i++)
-          os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
+    if (false)
+      for (int i = 0; i < rPE.nOutlinePoints(); i++)
+        os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
   }
 }
 
@@ -319,11 +319,11 @@ Phantom::printDefinitions (std::ostringstream& os) const
 
 
 /* NAME
- *   show              Show vector outline of Phantom to user
- *
- * SYNOPSIS
- *   show (pic)
- */
+*   show               Show vector outline of Phantom to user
+*
+* SYNOPSIS
+*   show (pic)
+*/
 
 #ifdef HAVE_SGP
 void 
@@ -331,9 +331,9 @@ Phantom::show () const
 {
   SGPDriver driverSGP ("Phantom Show");
   SGP sgp (driverSGP);
-
+  
   show (sgp);
-
+  
   std::cout << "Press return to continue";
   cio_kb_getc();
 }
@@ -343,26 +343,26 @@ Phantom::show (SGP& sgp) const
 {
   double wsize = m_xmax - m_xmin;
   if ((m_ymax - m_ymin) > wsize) 
-      wsize = m_ymax - m_ymin;
+    wsize = m_ymax - m_ymin;
   wsize *= 1.01;
   double halfWindow = wsize / 2;
-
+  
   double xcent = m_xmin + (m_xmax - m_xmin) / 2;
   double ycent = m_ymin + (m_ymax - m_ymin) / 2;
-
+  
   sgp.setWindow (xcent - halfWindow, ycent - halfWindow, xcent + halfWindow, ycent + halfWindow);
-
+  
   draw (sgp);
 }
 #endif
 
 
 /* NAME
- *   draw              Draw vector outline of Phantom
- *
- * SYNOPSIS
- *   draw ()
- */
+*   draw               Draw vector outline of Phantom
+*
+* SYNOPSIS
+*   draw ()
+*/
 
 #ifdef HAVE_SGP
 void 
@@ -375,13 +375,13 @@ Phantom::draw (SGP& sgp) const
 
 
 /* NAME
- *   addStdSheppLogan  Make head phantom of Shepp-Logan
- *
- * REFERENCES
- *   S. W. Rowland, "Computer Implementation of Image Reconstruction
- *     Formulas", in "Image Reconstruction from Projections: Implementation
- *     and Applications", edited by G. T. Herman, 1978.
- */
+*   addStdSheppLogan   Make head phantom of Shepp-Logan
+*
+* REFERENCES
+*   S. W. Rowland, "Computer Implementation of Image Reconstruction
+     Formulas", in "Image Reconstruction from Projections: Implementation
+     and Applications", edited by G. T. Herman, 1978.
+*/
 
 void 
 Phantom::addStdSheppLogan ()
@@ -401,12 +401,12 @@ Phantom::addStdSheppLogan ()
 
 
 /* NAME
- *   addStdHerman                      Standard head phantom of G. T. Herman
- *
- * REFERENCES
- *   G. T. Herman, "Image Reconstructions from Projections:  The Fundementals
- *     of Computed Tomography", 1979.
- */
+*   addStdHerman                       Standard head phantom of G. T. Herman
+*
+* REFERENCES
+*   G. T. Herman, "Image Reconstructions from Projections:  The Fundementals
+     of Computed Tomography", 1979.
+*/
 
 void 
 Phantom::addStdHerman ()
@@ -431,15 +431,15 @@ Phantom::addStdHerman ()
 
 
 /* NAME
- *    convertToImagefile               Make image array from Phantom
- *
- * SYNOPSIS
- *    pic_to_imagefile (pic, im, nsample)
- *    Phantom& pic             Phantom definitions
- *    ImageFile  *im           Computed pixel array
- *    int nsample              Number of samples along each axis for each pixel
- *                             (total samples per pixel = nsample * nsample)
- */
+*    convertToImagefile                Make image array from Phantom
+*
+* SYNOPSIS
+*    pic_to_imagefile (pic, im, nsample)
+*    Phantom& pic              Phantom definitions
+*    ImageFile  *im            Computed pixel array
+*    int nsample               Number of samples along each axis for each pixel
+                             (total samples per pixel = nsample * nsample)
+*/
 
 void
 Phantom::convertToImagefile (ImageFile& im, double dViewRatio, const int in_nsample, const int trace) const
@@ -448,84 +448,92 @@ Phantom::convertToImagefile (ImageFile& im, double dViewRatio, const int in_nsam
 }
 
 void 
-Phantom::convertToImagefile (ImageFile& im, const double dViewRatio, const int in_nsample, const int trace, const int colStart, const int colCount, bool bStoreAtColumnPos) const
+Phantom::convertToImagefile (ImageFile& im, const double dViewRatio, const int in_nsample, const int trace, 
+                             const int colStart, const int colCount, bool bStoreAtColumnPos) const
 {
-  int nx = im.nx();
-  int ny = im.ny();
-  if (nx < 2 || ny < 2)
-      return;
+  int iStorageOffset = (bStoreAtColumnPos ? colStart : 0);
+  convertToImagefile (im, im.nx(), dViewRatio, in_nsample, trace, colStart, colCount, iStorageOffset);
+}
 
+void 
+Phantom::convertToImagefile (ImageFile& im, const int iTotalRasterCols, const double dViewRatio, 
+            const int in_nsample, const int trace, const int colStart, const int colCount, int iStorageOffset) const
+{
+  const int nx = im.nx();
+  const int ny = im.ny();
+  if (nx < 2 || ny < 2)
+    return;
+  
   int nsample = in_nsample;
   if (nsample < 1)  
     nsample = 1;
-
+  
   double dx = m_xmax - m_xmin;
   double dy = m_ymax - m_ymin;
   double xcent = m_xmin + dx / 2;
   double ycent = m_ymin + dy / 2;
   double dHalflen = dViewRatio * (getDiameterBoundaryCircle() / SQRT2 / 2);
-
+  
   double xmin = xcent - dHalflen;
   double xmax = xcent + dHalflen;
   double ymin = ycent - dHalflen;
   double ymax = ycent + dHalflen;
-
+  
   // Each pixel holds the average of the intensity of the cell with (ix,iy) at the center of the pixel
   // Set major increments so that the last cell v[nx-1][ny-1] will start at xmax - xinc, ymax - yinc).
   // Set minor increments so that sample points are centered in cell
-
-  double xinc = (xmax - xmin) / nx;
+  
+  double xinc = (xmax - xmin) / (iTotalRasterCols);
   double yinc = (ymax - ymin) / ny;
-
+  
   double kxinc = xinc / nsample;               /* interval between samples */
   double kyinc = yinc / nsample;
   double kxofs = kxinc / 2;            /* offset of 1st point */
   double kyofs = kyinc / 2;
-
+  
   im.setAxisExtent (xmin, xmax, ymin, ymax);
   im.setAxisIncrement (xinc, yinc);
-
+  
   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;
-    }
-
+  
+  for (int ix = 0; ix < colCount; ix++) {
+    int iColStore = ix + iStorageOffset;
+    ImageFileColumn vCol = v[iColStore];
+    for (int iy = 0; iy < ny; iy++)
+      *vCol++ = 0;
+  }
+  
   double x_start = xmin + (colStart * xinc);
   for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
-      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 */
+    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 + iStorageOffset;
+      ImageFileColumn vCol = v[iColStore];
+      for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
+        double dAtten = 0;
+        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))
+              dAtten += rPElem.atten();
+        } // for kx
+        *vCol++ += dAtten;
+      } /* for iy */
+    }  /* for ix */
   }  /* for pelem */
   
-
+  
   if (nsample > 1) {
     double factor = 1.0 / static_cast<double>(nsample * nsample);
-
-
+    
+    
     for (int ix = 0; ix < colCount; ix++) {
-               int iColStore = ix;
-               if (bStoreAtColumnPos)
-                       iColStore += colStart;
-               for (int iy = 0; iy < ny; iy++)
-                       v[iColStore][iy] *= factor;
-       }
+      int iColStore = ix + iStorageOffset;
+      ImageFileColumn vCol = v[iColStore];
+      for (int iy = 0; iy < ny; iy++)
+        *vCol++ *= factor;
+    }
   }
 }
 
@@ -540,15 +548,15 @@ Phantom::convertToImagefile (ImageFile& im, const double dViewRatio, const int i
 
 
 PhantomElement::PhantomElement (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
-    : m_cx(cx), m_cy(cy), m_u(u), m_v(v), m_atten(atten), m_nPoints(0), m_xOutline(0), m_yOutline(0)
+: m_cx(cx), m_cy(cy), m_u(u), m_v(v), m_atten(atten), m_nPoints(0), m_xOutline(0), m_yOutline(0)
 {
   m_rot = convertDegreesToRadians (rot);   // convert angle to radians
-
+  
   m_type = convertNameToType (type);
-
+  
   makeTransformMatrices ();     // calc transform matrices between phantom and normalized phantomelement
   makeVectorOutline ();                // calculate vector outline of pelem 
-
+  
   m_rectLimits[0] = m_xmin;   m_rectLimits[1] = m_ymin;
   m_rectLimits[2] = m_xmax;   m_rectLimits[3] = m_ymax;
 }
@@ -557,8 +565,8 @@ PhantomElement::PhantomElement (const char *type, const double cx, const double
 
 PhantomElement::~PhantomElement ()
 {
-    delete m_xOutline;
-    delete m_yOutline;
+  delete m_xOutline;
+  delete m_yOutline;
 }
 
 void
@@ -578,29 +586,29 @@ PhantomElement::printDefinition (std::ostringstream& os) const
 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);
+  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);
 }
 
 const char* const
 PhantomElement::convertTypeToName (PhmElemType iType)
 {
   static char* pszType = "Unknown";
-
+  
   if (iType == PELEM_RECTANGLE)
     pszType = "rectangle";
   else if (iType == PELEM_TRIANGLE)
@@ -611,7 +619,7 @@ PhantomElement::convertTypeToName (PhmElemType iType)
     pszType = "sector";
   else if (iType == PELEM_SEGMENT)
     pszType = "segment";
-
+  
   return pszType;
 }
 
@@ -620,23 +628,23 @@ void
 PhantomElement::makeTransformMatrices ()
 {
   GRFMTX_2D temp;
-
+  
   // To map normalized Pelem coords to world Phantom 
   //     scale by (u, v)                                      
   //     rotate by rot                                 
   //     translate by (cx, cy)                        
-
+  
   scale_mtx2 (m_xformObjToPhm, m_u, m_v);
   rot_mtx2  (temp, m_rot);
   mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
   xlat_mtx2 (temp, m_cx, m_cy);
   mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
-
+  
   // to map world Phantom coodinates to normalized PElem coords
   //     translate by (-cx, -cy)
   //     rotate by -rot
   //     scale by (1/u, 1/v)
-
+  
   xlat_mtx2 (m_xformPhmToObj, -m_cx, -m_cy);
   rot_mtx2  (temp, -m_rot);
   mult_mtx2 (m_xformPhmToObj, temp, m_xformPhmToObj);
@@ -646,15 +654,15 @@ PhantomElement::makeTransformMatrices ()
 
 
 /* NAME
- *   pelem_make_points         INTERNAL routine to calculate point array for an pelem
- *
- * SYNOPSIS
- *   makepelempts (pelem)
- *   PELEM *pelem      pelem whose points we are calculating
- *
- * NOTES
- *   Called by phm_add_pelem()
- */
+*   pelem_make_points          INTERNAL routine to calculate point array for an pelem
+*
+* SYNOPSIS
+*   makepelempts (pelem)
+*   PELEM *pelem       pelem whose points we are calculating
+*
+* NOTES
+*   Called by phm_add_pelem()
+*/
 
 void
 PhantomElement::makeVectorOutline ()
@@ -662,7 +670,7 @@ PhantomElement::makeVectorOutline ()
   double radius, theta, start, stop;
   double xfact, yfact;
   int cpts;
-
+  
   m_nPoints = 0;
   switch (m_type) {
   case PELEM_RECTANGLE:
@@ -735,10 +743,10 @@ PhantomElement::makeVectorOutline ()
   
   // increase pelem extent by SCALE_PELEM_EXTENT to eliminate chance of
   //   missing actual pelem maximum due to polygonal sampling 
-
+  
   xfact = (m_xmax - m_xmin) * SCALE_PELEM_EXTENT;
   yfact = (m_ymax - m_ymin) * SCALE_PELEM_EXTENT;
-
+  
   m_xmin -= xfact;
   m_ymin -= yfact;
   m_xmax += xfact;
@@ -747,39 +755,39 @@ PhantomElement::makeVectorOutline ()
 
 
 /* NAME
- *   calc_arc                  Calculate outline of a arc of a circle
- *
- * SYNOPSIS
- *   calc_arc (x, y, xcent, ycent, pts, r, start, stop)
- *   double x[], y[];          Array of points
- *   int pts                   Number of points in array
- *   double xcent, ycent       Center of cirlce
- *   double r                  Radius of circle
- *   double start, stop                Beginning & ending angles
- */
+*   calc_arc                   Calculate outline of a arc of a circle
+*
+* SYNOPSIS
+*   calc_arc (x, y, xcent, ycent, pts, r, start, stop)
+*   double x[], y[];           Array of points
+*   int pts                    Number of points in array
+*   double xcent, ycent        Center of cirlce
+*   double r                   Radius of circle
+*   double start, stop         Beginning & ending angles
+*/
 
 void 
 PhantomElement::calcArcPoints (double x[], double y[], const int pts, const double xcent, const double ycent, const double r, const double start, const double stop)
 {
-    if (r <= 0.0)
-       sys_error (ERR_WARNING, "negative or zero radius in calc_arc()");
-
-    double theta = (stop - start) / (pts - 1); // angle incr. between points 
-    double c = cos(theta);
-    double s = sin(theta);
+  if (r <= 0.0)
+    sys_error (ERR_WARNING, "negative or zero radius in calc_arc()");
   
-    x[0] = r * cos (start) + xcent;
-    y[0] = r * sin (start) + ycent;
-
-    double xp = x[0] - xcent;
-    double yp = y[0] - ycent;
-    for (int i = 1; i < pts; i++) {
-       double xc = c * xp - s * yp;
-       double yc = s * xp + c * yp;
-       x[i] = xc + xcent;
-       y[i] = yc + ycent;
-       xp = xc;  yp = yc;
-    }
+  double theta = (stop - start) / (pts - 1);   // angle incr. between points 
+  double c = cos(theta);
+  double s = sin(theta);
+  
+  x[0] = r * cos (start) + xcent;
+  y[0] = r * sin (start) + ycent;
+  
+  double xp = x[0] - xcent;
+  double yp = y[0] - ycent;
+  for (int i = 1; i < pts; i++) {
+    double xc = c * xp - s * yp;
+    double yc = s * xp + c * yp;
+    x[i] = xc + xcent;
+    y[i] = yc + ycent;
+    xp = xc;  yp = yc;
+  }
 }
 
 
@@ -794,26 +802,26 @@ PhantomElement::calcArcPoints (double x[], double y[], const int pts, const doub
 void 
 PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const double u, const double v)
 {
-    calcArcPoints (x, y, m_nPoints, 0.0, 0.0, 1.0, 0.0, TWOPI);   // make a unit circle 
-    scale2d (x, y, m_nPoints, m_u, m_v);                            // scale to ellipse 
+  calcArcPoints (x, y, m_nPoints, 0.0, 0.0, 1.0, 0.0, TWOPI);   // make a unit circle 
+  scale2d (x, y, m_nPoints, m_u, m_v);                      // scale to ellipse 
 }
 
 
 /* NAME
- *   circle_pts                Calculate number of points to use for circle segment
- *
- * SYNOPSIS
- *   n = circle_pts (theta)
- *   int n             Number of points to use for arc
- *   double theta      Length of arc in radians
- */
+*   circle_pts         Calculate number of points to use for circle segment
+*
+* SYNOPSIS
+*   n = circle_pts (theta)
+*   int n              Number of points to use for arc
+*   double theta       Length of arc in radians
+*/
 
 int 
 PhantomElement::numCirclePoints (double theta)
 {
-    theta = clamp (theta, 0., TWOPI);
-
-    return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
+  theta = clamp (theta, 0., TWOPI);
+  
+  return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
 }
 
 
@@ -824,42 +832,42 @@ PhantomElement::clipLineWorldCoords (double& x1, double& y1, double& x2, double
   double cx1 = x1, cy1 = y1, cx2 = x2, cy2 = y2;
   if (! clip_rect (cx1, cy1, cx2, cy2, m_rectLimits))
     return false;
-       
+  
   // convert phantom coordinates to pelem coordinates 
   xform_mtx2 (m_xformPhmToObj, x1, y1);
   xform_mtx2 (m_xformPhmToObj, x2, y2);
-       
+  
   if (! clipLineNormalizedCoords (x1, y1, x2, y2))
     return false;
-
+  
   // convert standard pelem coordinates back to phantom coordinates 
   xform_mtx2 (m_xformObjToPhm, x1, y1);
   xform_mtx2 (m_xformObjToPhm, x2, y2);
-
+  
   return true;
 }
 
 
 /* NAME
- *   pelem_clip_line                   Clip pelem against an arbitrary line
- *
- * SYNOPSIS
- *   pelem_clip_line (pelem, x1, y1, x2, y2)
- *   PhantomElement& pelem;            Pelem to be clipped
- *   double *x1, *y1, *x2, *y2 Endpoints of line to be clipped
- *
- * RETURNS
- *   true   if line passes through pelem
- *             (x1, y1, x2, y2 hold coordinates of new line)
- *   false  if line do not pass through pelem
- *             (x1, y1, x2, y2 are undefined)
- */
+*   pelem_clip_line                    Clip pelem against an arbitrary line
+*
+* SYNOPSIS
+*   pelem_clip_line (pelem, x1, y1, x2, y2)
+*   PhantomElement& pelem;             Pelem to be clipped
+*   double *x1, *y1, *x2, *y2  Endpoints of line to be clipped
+*
+* RETURNS
+*   true   if line passes through pelem
+             (x1, y1, x2, y2 hold coordinates of new line)
+*   false  if line do not pass through pelem
+             (x1, y1, x2, y2 are undefined)
+*/
 
 bool
 PhantomElement::clipLineNormalizedCoords (double& x1, double& y1, double& x2, double& y2) const
 {
   bool accept = false;
-
+  
   switch (m_type) {
   case PELEM_RECTANGLE:
     double rect[4];
@@ -883,7 +891,7 @@ PhantomElement::clipLineNormalizedCoords (double& x1, double& y1, double& x2, do
     sys_error (ERR_WARNING, "Illegal pelem type %d [pelem_clip_line]", m_type);
     break;
   }
-
+  
   return(accept);
 }
 
@@ -909,7 +917,7 @@ PhantomElement::isPointInside (double x, double y, const CoordType coord_type) c
     sys_error(ERR_WARNING, "Illegal coordinate type in pelem_is_point_inside");
     return (false);
   }
-
+  
   switch (m_type) {
   case PELEM_RECTANGLE:
     if (x > 1. || x < -1. || y > 1. || y < -1.)
@@ -931,13 +939,13 @@ PhantomElement::isPointInside (double x, double y, const CoordType coord_type) c
     else
       return (true);
     break;
-
+    
     // for clipping segments & sectors, must NOT scale by (1/u, 1/v)
     // because this destroys information about size of arc component 
-
+    
   case PELEM_SEGMENT:
     if (x > 1. || x < -1. || y > 0.)
-       return (false);         // clip against y > 0 
+      return (false);          // clip against y > 0 
     x *= m_u;                  // put back u & v scale 
     y *= m_v;
     if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
@@ -946,22 +954,22 @@ PhantomElement::isPointInside (double x, double y, const CoordType coord_type) c
       return (true);
     break;
   case PELEM_SECTOR:
-      if (x > 1. || x < -1. || y > 1.)   // extent 
+    if (x > 1. || x < -1. || y > 1.)   // extent 
       return (false);
-      if (y > 1. - x || y > 1. + x)      // triangle   
-         return (false);                      // clip against triangle 
-      x *= m_u;                       // circle: put back u & v scale 
+    if (y > 1. - x || y > 1. + x)      // triangle     
+      return (false);                 // clip against triangle 
+    x *= m_u;                 // circle: put back u & v scale 
     y *= m_v;
     if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
-       return (false);                // clip against circle 
+      return (false);                 // clip against circle 
     else
       return (true);
-  break;
+    break;
   default:
     sys_error (ERR_WARNING, "Illegal pelem type in pelem_is_point_inside()");
     break;
   }
-
+  
   return (false);
 }