libctsim/phantom.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **     Name:                   phm.cpp
   5 **     Purpose:                Routines for phantom objects
   6 **     Progammer:              Kevin Rosenberg
   7 **     Date Started:           Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: phantom.cpp,v 1.26 2001/01/02 16:02:13 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30 const int PhantomElement::POINTS_PER_CIRCLE = 360;
  31 const double PhantomElement::SCALE_PELEM_EXTENT=0.005;  // increase pelem limits by 0.5%
  32
  33 const int Phantom::PHM_INVALID = -1;
  34 const int Phantom::PHM_HERMAN = 0;
  35 const int Phantom::PHM_B_HERMAN = 1;
  36 const int Phantom::PHM_SHEPP_LOGAN = 2;
  37 const int Phantom::PHM_B_SHEPP_LOGAN = 3;
  38 const int Phantom::PHM_UNITPULSE = 4;
  39
  40 const char* Phantom::s_aszPhantomName[] =
  41 {
  42   {"herman"},
  43   {"herman-b"},
  44   {"shepp-logan"},
  45   {"shepp-logan-b"},
  46   {"unitpulse"},
  47 };
  48
  49 const char* Phantom::s_aszPhantomTitle[] =
  50 {
  51   {"Herman Head"},
  52   {"Herman Head (Bordered)"},
  53   {"Shepp-Logan"},
  54   {"Shepp-Logan (Bordered)"},
  55   {"Unit Pulse"},
  56 };
  57
  58 const int Phantom::s_iPhantomCount = sizeof(s_aszPhantomName) / sizeof(const char*);
  59
  60
  61 // CLASS IDENTIFICATION
  62 //   Phantom
  63 //
  64
  65 Phantom::Phantom ()
  66 {
  67   init ();
  68 }
  69
  70
  71 Phantom::Phantom (const char* const phmName)
  72 {
  73   init ();
  74   createFromPhantom (phmName);
  75 }
  76
  77 void
  78 Phantom::init ()
  79 {
  80   m_nPElem = 0;
  81   m_xmin = 1E30;
  82   m_xmax = -1E30;
  83   m_ymin = 1E30;
  84   m_ymax = -1E30;
  85   m_diameter = 0;
  86   m_composition = P_PELEMS;
  87   m_fail = false;
  88   m_id = PHM_INVALID;
  89 }
  90
  91 Phantom::~Phantom ()
  92 {
  93   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
  94     delete *i;
  95   }
  96 }
  97
  98
  99 const char*
 100 Phantom::convertPhantomIDToName (int phmID)
 101 {
 102   static const char *name = "";
 103
 104   if (phmID >= 0 && phmID < s_iPhantomCount)
 105       return (s_aszPhantomName[phmID]);
 106
 107   return (name);
 108 }
 109
 110 const char*
 111 Phantom::convertPhantomIDToTitle (int phmID)
 112 {
 113   static const char *title = "";
 114
 115   if (phmID >= 0 && phmID < s_iPhantomCount)
 116       return (s_aszPhantomName[phmID]);
 117
 118   return (title);
 119 }
 120
 121 int
 122 Phantom::convertNameToPhantomID (const char* const phmName)
 123 {
 124   int id = PHM_INVALID;
 125
 126   for (int i = 0; i < s_iPhantomCount; i++)
 127       if (strcasecmp (phmName, s_aszPhantomName[i]) == 0) {
 128           id = i;
 129           break;
 130       }
 131
 132   return (id);
 133 }
 134
 135
 136 bool
 137 Phantom::createFromPhantom (const char* const phmName)
 138 {
 139   int phmid = convertNameToPhantomID (phmName);
 140   if (phmid == PHM_INVALID) {
 141     m_fail = true;
 142     m_failMessage = "Invalid phantom name ";
 143     m_failMessage += phmName;
 144     return false;
 145   }
 146
 147   m_name = phmName;
 148   createFromPhantom (phmid);
 149   return true;
 150 }
 151
 152 bool
 153 Phantom::createFromPhantom (const int phmid)
 154 {
 155   switch (phmid)
 156     {
 157     case PHM_HERMAN:
 158       addStdHerman();
 159       break;
 160     case PHM_B_HERMAN:
 161       addStdHermanBordered();
 162       break;
 163     case PHM_SHEPP_LOGAN:
 164       addStdSheppLogan();
 165       break;
 166     case PHM_B_SHEPP_LOGAN:
 167       addStdSheppLoganBordered();
 168       break;
 169     case PHM_UNITPULSE:
 170       m_composition = P_UNIT_PULSE;
 171       addPElem ("rectangle", 0., 0., 100., 100., 0., 0.);     // outline
 172       addPElem ("ellipse", 0., 0., 1., 1., 0., 1.);           // pulse
 173       break;
 174     default:
 175       m_fail = true;
 176       m_failMessage = "Illegal phantom id ";
 177       m_failMessage += phmid;
 178       return false;
 179     }
 180
 181   m_id = phmid;
 182
 183   return true;
 184 }
 185
 186
 187 /* METHOD IDENTIFICATION
 188  *   createFromFile          Add PhantomElements from file
 189  *
 190  * SYNOPSIS
 191  *   createFromFile (filename)
 192  *
 193  * RETURNS
 194  *   true if pelem were added
 195  *   false if an pelem not added
 196  */
 197
 198 bool
 199 Phantom::createFromFile (const char* const fname)
 200 {
 201   bool bGoodFile = true;
 202   FILE *fp;
 203
 204   if ((fp = fopen (fname, "r")) == NULL)
 205     return (false);
 206
 207   m_name = fname;
 208
 209   while (1) {
 210     double cx, cy, u, v, rot, dens;
 211     char pelemtype[80];
 212
 213     int status = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf", pelemtype, &cx, &cy, &u, &v, &rot, &dens);
 214
 215     if (status == static_cast<int>(EOF))
 216       break;
 217     else if (status != 7) {
 218       sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromFile]", fname);
 219       bGoodFile = false;
 220     }
 221     addPElem (pelemtype, cx, cy, u, v, rot, dens);
 222   }
 223
 224   fclose (fp);
 225
 226   return (bGoodFile);
 227 }
 228
 229 bool
 230 Phantom::fileWrite (const char* const fname)
 231 {
 232   fstream file (fname, ios::out);
 233
 234   if (! file.fail())
 235     printDefinitions (file);
 236   return ! file.fail();
 237 }
 238
 239 /* NAME
 240  *   addPElem           Add pelem
 241  *
 242  * SYNOPSIS
 243  *   addPElem (type, cx, cy, u, v, rot, atten)
 244  *   char *type         type of pelem (box, ellipse, etc)
 245  *   double cx, cy      pelem center
 246  *   double u,v         pelem size
 247  *   double rot         rotation angle of pelem (in degrees)
 248  *   double atten       x-ray attenuation cooefficient
 249  */
 250
 251 void
 252 Phantom::addPElem (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
 253 {
 254   PhantomElement *pelem = new PhantomElement (type, cx, cy, u, v, rot, atten);
 255
 256   m_listPElem.push_front (pelem);
 257
 258   // update phantom limits
 259   if (m_xmin > pelem->xmin())    m_xmin = pelem->xmin();
 260   if (m_xmax < pelem->xmax())    m_xmax = pelem->xmax();
 261   if (m_ymin > pelem->ymin())    m_ymin = pelem->ymin();
 262   if (m_ymax < pelem->ymax())    m_ymax = pelem->ymax();
 263
 264   if (m_diameter < pelem->diameter())
 265     m_diameter = pelem->diameter();
 266
 267   //  m_diameter = lineLength(m_xmin, m_ymin, m_xmax, m_ymax);
 268
 269   m_nPElem++;
 270 }
 271
 272
 273 /*----------------------------------------------------------------------*/
 274 /*                      Input-Output Routines                           */
 275 /*----------------------------------------------------------------------*/
 276
 277
 278 /* NAME
 279  *   print                              Print vertices of Phantom pelems
 280  *
 281  * SYNOPSIS
 282  *   print (phm)
 283  */
 284
 285 void
 286 Phantom::print (std::ostream& os) const
 287 {
 288   os << "Number of PElements: " << m_nPElem << "\n";
 289   os << "Limits: xmin=" << m_xmin << ", ymin=" << m_ymin << ", xmax=" << m_xmax << ", ymax=" << m_ymax << "\n";
 290
 291   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
 292     const PhantomElement& rPE = **i;
 293       os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
 294       os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
 295       os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
 296
 297       if (false)
 298         for (int i = 0; i < rPE.nOutlinePoints(); i++)
 299           os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
 300   }
 301 }
 302 void
 303 Phantom::print (std::ostringstream& os) const
 304 {
 305   os << "Number of PElements: " << m_nPElem << "\n";
 306   os << "Limits: xmin=" << m_xmin << ", ymin=" << m_ymin << ", xmax=" << m_xmax << ", ymax=" << m_ymax << "\n";
 307
 308   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
 309     const PhantomElement& rPE = **i;
 310       os << "PhantomElement: nPoints=" << rPE.nOutlinePoints();
 311       os << ", atten=" << rPE.atten() << " rot=" << convertRadiansToDegrees (rPE.rot()) << "\n";
 312       os << "xmin=" << rPE.xmin() << ", ymin=" << rPE.ymin() << ", xmax=" << rPE.xmax() << ", ymax=" << rPE.ymax() << "\n";
 313
 314       if (false)
 315         for (int i = 0; i < rPE.nOutlinePoints(); i++)
 316           os << rPE.xOutline()[i] << "," << rPE.yOutline()[i] << "\n";
 317   }
 318 }
 319
 320 void
 321 Phantom::printDefinitions (std::ostream& os) const
 322 {
 323   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
 324     const PhantomElement& rPE = **i;
 325     rPE.printDefinition (os);
 326   }
 327 }
 328
 329 void
 330 Phantom::printDefinitions (std::ostringstream& os) const
 331 {
 332   for (PElemConstIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
 333     const PhantomElement& rPE = **i;
 334     rPE.printDefinition (os);
 335   }
 336 }
 337
 338
 339 /* NAME
 340  *   show               Show vector outline of Phantom to user
 341  *
 342  * SYNOPSIS
 343  *   show (pic)
 344  */
 345
 346 #ifdef HAVE_SGP
 347 void
 348 Phantom::show () const
 349 {
 350   SGPDriver driverSGP ("Phantom Show");
 351   SGP sgp (driverSGP);
 352
 353   show (sgp);
 354
 355   std::cout << "Press return to continue";
 356   cio_kb_getc();
 357 }
 358
 359 void
 360 Phantom::show (SGP& sgp) const
 361 {
 362   double wsize = m_xmax - m_xmin;
 363   if ((m_ymax - m_ymin) > wsize)
 364       wsize = m_ymax - m_ymin;
 365   wsize *= 1.01;
 366   double halfWindow = wsize / 2;
 367
 368   double xcent = m_xmin + (m_xmax - m_xmin) / 2;
 369   double ycent = m_ymin + (m_ymax - m_ymin) / 2;
 370
 371   sgp.setWindow (xcent - halfWindow, ycent - halfWindow, xcent + halfWindow, ycent + halfWindow);
 372
 373   draw (sgp);
 374 }
 375 #endif
 376
 377
 378 /* NAME
 379  *   draw               Draw vector outline of Phantom
 380  *
 381  * SYNOPSIS
 382  *   draw ()
 383  */
 384
 385 #ifdef HAVE_SGP
 386 void
 387 Phantom::draw (SGP& sgp) const
 388 {
 389   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
 390     sgp.polylineAbs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
 391 }
 392 #endif
 393
 394
 395 /* NAME
 396  *   addStdSheppLogan   Make head phantom of Shepp-Logan
 397  *
 398  * REFERENCES
 399  *   S. W. Rowland, "Computer Implementation of Image Reconstruction
 400  *      Formulas", in "Image Reconstruction from Projections: Implementation
 401  *      and Applications", edited by G. T. Herman, 1978.
 402  */
 403
 404 void
 405 Phantom::addStdSheppLogan ()
 406 {
 407   addPElem ("ellipse",  0.0000,  0.0000, 0.6900,  0.9200,   0.0,  1.00);
 408   addPElem ("ellipse",  0.0000, -0.0184, 0.6624,  0.8740,   0.0, -0.98);
 409   addPElem ("ellipse",  0.2200,  0.0000, 0.1100,  0.3100, -18.0, -0.02);
 410   addPElem ("ellipse", -0.2200,  0.0000, 0.1600,  0.4100,  18.0, -0.02);
 411   addPElem ("ellipse",  0.0000,  0.3500, 0.2100,  0.2500,   0.0,  0.01);
 412   addPElem ("ellipse",  0.0000,  0.1000, 0.0460,  0.0460,   0.0,  0.01);
 413   addPElem ("ellipse",  0.0000, -0.1000, 0.0460,  0.0460,   0.0,  0.01);
 414   addPElem ("ellipse", -0.0800, -0.6050, 0.0460,  0.0230,   0.0,  0.01);
 415   addPElem ("ellipse",  0.0000, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
 416   addPElem ("ellipse",  0.0600, -0.6050, 0.0230,  0.0230,   0.0,  0.01);
 417   addPElem ("ellipse",  0.5538, -0.3858, 0.0330,  0.2060, -18.0,  0.03);
 418 }
 419
 420 void
 421 Phantom::addStdSheppLoganBordered ()
 422 {
 423   addStdSheppLogan ();
 424   addPElem ("rectangle", 0.000, 0.0000, 0.8600, 1.150, 0.0, 0.00);
 425 }
 426
 427 /* NAME
 428  *   addStdHerman                       Standard head phantom of G. T. Herman
 429  *
 430  * REFERENCES
 431  *   G. T. Herman, "Image Reconstructions from Projections:  The Fundementals
 432  *      of Computed Tomography", 1979.
 433  */
 434
 435 void
 436 Phantom::addStdHerman ()
 437 {
 438   addPElem ("ellipse",  0.000,  1.50,  0.375, 0.3000,  90.00, -0.003);
 439   addPElem ("ellipse",  0.675, -0.75,  0.225, 0.1500, 140.00,  0.010);
 440   addPElem ("ellipse",  0.750,  1.50,  0.375, 0.2250,  50.00,  0.003);
 441   addPElem ("segment",  1.375, -7.50,  1.100, 0.6250,  19.20, -0.204);
 442   addPElem ("segment",  1.375, -7.50,  1.100, 4.3200,  19.21,  0.204);
 443   addPElem ("segment",  0.000, -2.25,  1.125, 0.3750,   0.00, -0.003);
 444   addPElem ("segment",  0.000, -2.25,  1.125, 3.0000,   0.00,  0.003);
 445   addPElem ("segment", -1.000,  3.75,  1.000, 0.5000, 135.00, -0.003);
 446   addPElem ("segment", -1.000,  3.75,  1.000, 3.0000, 135.00,  0.003);
 447   addPElem ("segment",  1.000,  3.75,  1.000, 0.5000, 225.00, -0.003);
 448   addPElem ("segment",  1.000,  3.75,  1.000, 3.0000, 225.00,  0.003);
 449   addPElem ("triangle", 5.025,  3.75,  1.125, 0.5000, 110.75,  0.206);
 450   addPElem ("triangle",-5.025,  3.75,  1.125, 0.9000,-110.75,  0.206);
 451   addPElem ("ellipse",  0.000,  0.00,  8.625, 6.4687,  90.00,  0.416);
 452   addPElem ("ellipse",  0.000,  0.00,  7.875, 5.7187,  90.00, -0.206);
 453 }
 454
 455 void
 456 Phantom::addStdHermanBordered ()
 457 {
 458   addStdHerman();
 459   addPElem ("rectangle",  0.000, 0.00, 10.780, 8.110,  90.00, 0.000);
 460 }
 461
 462
 463 /* NAME
 464  *    convertToImagefile                Make image array from Phantom
 465  *
 466  * SYNOPSIS
 467  *    pic_to_imagefile (pic, im, nsample)
 468  *    Phantom& pic              Phantom definitions
 469  *    ImageFile  *im            Computed pixel array
 470  *    int nsample               Number of samples along each axis for each pixel
 471  *                              (total samples per pixel = nsample * nsample)
 472  */
 473
 474 void
 475 Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const
 476 {
 477   convertToImagefile (im, in_nsample, trace, 0, im.nx(), true);
 478 }
 479
 480 void
 481 Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount, bool bStoreAtColumnPos) const
 482 {
 483   int nx = im.nx();
 484   int ny = im.ny();
 485   if (nx < 2 || ny < 2)
 486       return;
 487
 488   int nsample = in_nsample;
 489   if (nsample < 1)
 490     nsample = 1;
 491
 492   double dx = m_xmax - m_xmin;
 493   double dy = m_ymax - m_ymin;
 494   double xcent = m_xmin + dx / 2;
 495   double ycent = m_ymin + dy / 2;
 496   double phmlen = (dx > dy ? dx : dy);
 497
 498   double phmradius = phmlen / 2;
 499
 500   double xmin = xcent - phmradius;
 501   double xmax = xcent + phmradius;
 502   double ymin = ycent - phmradius;
 503   double ymax = ycent + phmradius;
 504
 505   // Each pixel holds the average of the intensity of the cell with (ix,iy) at the center of the pixel
 506   // Set major increments so that the last cell v[nx-1][ny-1] will start at xmax - xinc, ymax - yinc).
 507   // Set minor increments so that sample points are centered in cell
 508
 509   double xinc = (xmax - xmin) / nx;
 510   double yinc = (ymax - ymin) / ny;
 511
 512   double kxinc = xinc / nsample;                /* interval between samples */
 513   double kyinc = yinc / nsample;
 514   double kxofs = kxinc / 2;             /* offset of 1st point */
 515   double kyofs = kyinc / 2;
 516
 517   im.setAxisExtent (xmin, xmax, ymin, ymax);
 518   im.setAxisIncrement (xinc, yinc);
 519
 520   ImageFileArray v = im.getArray();
 521
 522   for (int ix = 0; ix < colCount; ix++)
 523     for (int iy = 0; iy < ny; iy++) {
 524       int iColStore = ix;
 525       if (bStoreAtColumnPos)
 526                 iColStore += colStart;
 527       v[iColStore][iy] = 0;
 528     }
 529
 530   double x_start = xmin + (colStart * xinc);
 531   for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
 532       const PhantomElement& rPElem = **pelem;
 533       double x, y, xi, yi;
 534       int ix, iy, kx, ky;
 535       for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) {
 536         int iColStore = ix;
 537         if (bStoreAtColumnPos)
 538           iColStore += colStart;
 539         for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
 540           for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
 541             for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
 542               if (rPElem.isPointInside (xi, yi, PHM_COORD) == TRUE)
 543                         v[iColStore][iy] += rPElem.atten();
 544               } // for kx
 545           } /* for iy */
 546       }  /* for ix */
 547   }  /* for pelem */
 548
 549
 550   if (nsample > 1) {
 551     double factor = 1.0 / static_cast<double>(nsample * nsample);
 552
 553
 554     for (int ix = 0; ix < colCount; ix++) {
 555                 int iColStore = ix;
 556                 if (bStoreAtColumnPos)
 557                         iColStore += colStart;
 558                 for (int iy = 0; iy < ny; iy++)
 559                         v[iColStore][iy] *= factor;
 560         }
 561   }
 562 }
 563
 564 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 565 // CLASS IDENTIFICATION
 566 //
 567 //      PhantomElement
 568 //
 569 // PURPOSE
 570 //
 571 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 572
 573
 574 PhantomElement::PhantomElement (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
 575     : m_cx(cx), m_cy(cy), m_u(u), m_v(v), m_atten(atten), m_nPoints(0), m_xOutline(0), m_yOutline(0)
 576 {
 577   m_rot = convertDegreesToRadians (rot);   // convert angle to radians
 578
 579   m_type = convertNameToType (type);
 580
 581   makeTransformMatrices ();     // calc transform matrices between phantom and normalized phantomelement
 582   makeVectorOutline ();         // calculate vector outline of pelem
 583
 584   // Find maximum diameter of Object
 585   double r2Max = 0;
 586   for (int i = 0; i < m_nPoints; i++) {
 587     double r2 = (m_xOutline[i] * m_xOutline[i]) + (m_yOutline[i] * m_yOutline[i]);
 588     if (r2 > r2Max)
 589       r2Max = r2;
 590   }
 591   m_diameter = 2 * sqrt( r2Max );
 592
 593   m_rectLimits[0] = m_xmin;   m_rectLimits[1] = m_ymin;
 594   m_rectLimits[2] = m_xmax;   m_rectLimits[3] = m_ymax;
 595 }
 596
 597
 598
 599 PhantomElement::~PhantomElement ()
 600 {
 601     delete m_xOutline;
 602     delete m_yOutline;
 603 }
 604
 605 void
 606 PhantomElement::printDefinition (std::ostream& os) const
 607 {
 608   os << convertTypeToName (m_type) << " " << m_cx << " " << m_cy << " " << m_u << " "
 609     << m_v << " " << convertRadiansToDegrees (m_rot) << " " << m_atten << "\n";
 610 }
 611
 612 void
 613 PhantomElement::printDefinition (std::ostringstream& os) const
 614 {
 615   os << convertTypeToName (m_type) << " " << m_cx << " " << m_cy << " " << m_u << " "
 616     << m_v << " " << convertRadiansToDegrees (m_rot) << " " << m_atten << "\n";
 617 }
 618
 619 PhmElemType
 620 PhantomElement::convertNameToType (const char* const typeName)
 621 {
 622     PhmElemType type = PELEM_INVALID;
 623
 624     if (strcasecmp (typeName, "rectangle") == 0)
 625         type = PELEM_RECTANGLE;
 626     else if (strcasecmp (typeName, "triangle") == 0)
 627         type = PELEM_TRIANGLE;
 628     else if (strcasecmp (typeName, "ellipse") == 0)
 629         type = PELEM_ELLIPSE;
 630     else if (strcasecmp (typeName, "sector") == 0)
 631         type = PELEM_SECTOR;
 632     else if (strcasecmp (typeName, "segment") == 0)
 633       type = PELEM_SEGMENT;
 634     else
 635         sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
 636
 637     return (type);
 638 }
 639
 640 const char* const
 641 PhantomElement::convertTypeToName (PhmElemType iType)
 642 {
 643   static char* pszType = "Unknown";
 644
 645   if (iType == PELEM_RECTANGLE)
 646     pszType = "rectangle";
 647   else if (iType == PELEM_TRIANGLE)
 648     pszType = "triangle";
 649   else if (iType == PELEM_ELLIPSE)
 650     pszType = "ellipse";
 651   else if (iType == PELEM_SECTOR)
 652     pszType = "sector";
 653   else if (iType == PELEM_SEGMENT)
 654     pszType = "segment";
 655
 656   return pszType;
 657 }
 658
 659
 660 void
 661 PhantomElement::makeTransformMatrices ()
 662 {
 663   GRFMTX_2D temp;
 664
 665   // To map normalized Pelem coords to world Phantom
 666   //     scale by (u, v)
 667   //     rotate by rot
 668   //     translate by (cx, cy)
 669
 670   scale_mtx2 (m_xformObjToPhm, m_u, m_v);
 671   rot_mtx2  (temp, m_rot);
 672   mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
 673   xlat_mtx2 (temp, m_cx, m_cy);
 674   mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
 675
 676   // to map world Phantom coodinates to normalized PElem coords
 677   //     translate by (-cx, -cy)
 678   //     rotate by -rot
 679   //     scale by (1/u, 1/v)
 680
 681   xlat_mtx2 (m_xformPhmToObj, -m_cx, -m_cy);
 682   rot_mtx2  (temp, -m_rot);
 683   mult_mtx2 (m_xformPhmToObj, temp, m_xformPhmToObj);
 684   scale_mtx2 (temp, 1 / m_u, 1 / m_v);
 685   mult_mtx2 (m_xformPhmToObj, temp, m_xformPhmToObj);
 686 }
 687
 688
 689 /* NAME
 690  *   pelem_make_points          INTERNAL routine to calculate point array for an pelem
 691  *
 692  * SYNOPSIS
 693  *   makepelempts (pelem)
 694  *   PELEM *pelem       pelem whose points we are calculating
 695  *
 696  * NOTES
 697  *   Called by phm_add_pelem()
 698  */
 699
 700 void
 701 PhantomElement::makeVectorOutline ()
 702 {
 703   double radius, theta, start, stop;
 704   double xfact, yfact;
 705   int cpts;
 706
 707   m_nPoints = 0;
 708   switch (m_type) {
 709   case PELEM_RECTANGLE:
 710     m_nPoints = 5;
 711     m_xOutline = new double [m_nPoints];
 712     m_yOutline = new double [m_nPoints];
 713     m_xOutline[0] =-m_u;        m_yOutline[0] =-m_v;
 714     m_xOutline[1] = m_u;        m_yOutline[1] =-m_v;
 715     m_xOutline[2] = m_u;        m_yOutline[2] = m_v;
 716     m_xOutline[3] =-m_u;        m_yOutline[3] = m_v;
 717     m_xOutline[4] =-m_u;        m_yOutline[4] =-m_v;
 718     break;
 719   case PELEM_TRIANGLE:
 720     m_nPoints = 4;
 721     m_xOutline = new double [m_nPoints];
 722     m_yOutline = new double [m_nPoints];
 723     m_xOutline[0] =-m_u;        m_yOutline[0] = 0.0;
 724     m_xOutline[1] = m_u;        m_yOutline[1] = 0.0;
 725     m_xOutline[2] = 0.0;        m_yOutline[2] = m_v;
 726     m_xOutline[3] =-m_u;        m_yOutline[3] = 0.0;
 727     break;
 728   case PELEM_ELLIPSE:
 729     cpts = numCirclePoints (TWOPI);
 730     m_nPoints = cpts;
 731     m_xOutline = new double [m_nPoints];
 732     m_yOutline = new double [m_nPoints];
 733     calcEllipsePoints (m_xOutline, m_yOutline, cpts, m_u, m_v);
 734     break;
 735   case PELEM_SECTOR:
 736     radius = sqrt(m_u * m_u + m_v * m_v);
 737     theta = atan(m_u / m_v);            // angle with y-axis
 738     start = 3.0 * HALFPI - theta;
 739     stop  = 3.0 * HALFPI + theta;
 740     cpts = numCirclePoints (stop - start);
 741     m_nPoints = 3 + cpts;
 742     m_xOutline = new double [m_nPoints];
 743     m_yOutline = new double [m_nPoints];
 744
 745     m_xOutline[0] = 0.0;                m_yOutline[0] = m_v;
 746     m_xOutline[1] =-m_u;                m_yOutline[1] = 0.0;
 747     calcArcPoints (&m_xOutline[2], &m_yOutline[2], cpts, 0.0, m_v, radius, start, stop);
 748     m_xOutline[cpts + 2] = 0.0;
 749     m_yOutline[cpts + 2] = m_v;
 750     break;
 751   case PELEM_SEGMENT:
 752     radius = sqrt(m_u * m_u + m_v * m_v);
 753     theta = atan (m_u / m_v);           // angle with y-axis
 754     start = 3.0 * HALFPI - theta;
 755     stop  = 3.0 * HALFPI + theta;
 756
 757     cpts = numCirclePoints (stop - start);
 758     m_nPoints = cpts + 1;
 759     m_xOutline = new double [m_nPoints];
 760     m_yOutline = new double [m_nPoints];
 761
 762     calcArcPoints (m_xOutline, m_yOutline, cpts, 0.0, m_v, radius, start, stop);
 763     m_xOutline[cpts] = -m_u;
 764     m_yOutline[cpts] = 0.0;
 765     break;
 766   default:
 767     sys_error(ERR_WARNING, "illegal pelem type %d [makeVectorOutline]", m_type);
 768     return;
 769   }
 770
 771   rotate2d (m_xOutline, m_yOutline, m_nPoints, m_rot);
 772   xlat2d (m_xOutline, m_yOutline, m_nPoints, m_cx, m_cy);
 773
 774   minmax_array (m_xOutline, m_nPoints, m_xmin, m_xmax);
 775   minmax_array (m_yOutline, m_nPoints, m_ymin, m_ymax);
 776
 777   // increase pelem extent by SCALE_PELEM_EXTENT to eliminate chance of
 778   //   missing actual pelem maximum due to polygonal sampling
 779
 780   xfact = (m_xmax - m_xmin) * SCALE_PELEM_EXTENT;
 781   yfact = (m_ymax - m_ymin) * SCALE_PELEM_EXTENT;
 782
 783   m_xmin -= xfact;
 784   m_ymin -= yfact;
 785   m_xmax += xfact;
 786   m_ymax += yfact;
 787 }
 788
 789
 790 /* NAME
 791  *   calc_arc                   Calculate outline of a arc of a circle
 792  *
 793  * SYNOPSIS
 794  *   calc_arc (x, y, xcent, ycent, pts, r, start, stop)
 795  *   double x[], y[];           Array of points
 796  *   int pts                    Number of points in array
 797  *   double xcent, ycent        Center of cirlce
 798  *   double r                   Radius of circle
 799  *   double start, stop         Beginning & ending angles
 800  */
 801
 802 void
 803 PhantomElement::calcArcPoints (double x[], double y[], const int pts, const double xcent, const double ycent, const double r, const double start, const double stop)
 804 {
 805     if (r <= 0.0)
 806         sys_error (ERR_WARNING, "negative or zero radius in calc_arc()");
 807
 808     double theta = (stop - start) / (pts - 1);  // angle incr. between points
 809     double c = cos(theta);
 810     double s = sin(theta);
 811
 812     x[0] = r * cos (start) + xcent;
 813     y[0] = r * sin (start) + ycent;
 814
 815     double xp = x[0] - xcent;
 816     double yp = y[0] - ycent;
 817     for (int i = 1; i < pts; i++) {
 818         double xc = c * xp - s * yp;
 819         double yc = s * xp + c * yp;
 820         x[i] = xc + xcent;
 821         y[i] = yc + ycent;
 822         xp = xc;  yp = yc;
 823     }
 824 }
 825
 826
 827 // NAME
 828 //   PhantomElement::calcEllipsePoints    Calculate outline of a ellipse
 829 //
 830 // SYNOPSIS
 831 //   calcEllipsePoints ()
 832 //
 833
 834
 835 void
 836 PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const double u, const double v)
 837 {
 838     calcArcPoints (x, y, m_nPoints, 0.0, 0.0, 1.0, 0.0, TWOPI);   // make a unit circle
 839     scale2d (x, y, m_nPoints, m_u, m_v);                             // scale to ellipse
 840 }
 841
 842
 843 /* NAME
 844  *   circle_pts         Calculate number of points to use for circle segment
 845  *
 846  * SYNOPSIS
 847  *   n = circle_pts (theta)
 848  *   int n              Number of points to use for arc
 849  *   double theta       Length of arc in radians
 850  */
 851
 852 int
 853 PhantomElement::numCirclePoints (double theta)
 854 {
 855     theta = clamp (theta, 0., TWOPI);
 856
 857     return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
 858 }
 859
 860
 861 bool
 862 PhantomElement::clipLineWorldCoords (double& x1, double& y1, double& x2, double &y2) const
 863 {
 864   /* check if ray is outside of pelem extents */
 865   double cx1 = x1, cy1 = y1, cx2 = x2, cy2 = y2;
 866   if (! clip_rect (cx1, cy1, cx2, cy2, m_rectLimits))
 867     return false;
 868
 869   // convert phantom coordinates to pelem coordinates
 870   xform_mtx2 (m_xformPhmToObj, x1, y1);
 871   xform_mtx2 (m_xformPhmToObj, x2, y2);
 872
 873   if (! clipLineNormalizedCoords (x1, y1, x2, y2))
 874     return false;
 875
 876   // convert standard pelem coordinates back to phantom coordinates
 877   xform_mtx2 (m_xformObjToPhm, x1, y1);
 878   xform_mtx2 (m_xformObjToPhm, x2, y2);
 879
 880   return true;
 881 }
 882
 883
 884 /* NAME
 885  *   pelem_clip_line                    Clip pelem against an arbitrary line
 886  *
 887  * SYNOPSIS
 888  *   pelem_clip_line (pelem, x1, y1, x2, y2)
 889  *   PhantomElement& pelem;             Pelem to be clipped
 890  *   double *x1, *y1, *x2, *y2  Endpoints of line to be clipped
 891  *
 892  * RETURNS
 893  *   true   if line passes through pelem
 894  *              (x1, y1, x2, y2 hold coordinates of new line)
 895  *   false  if line do not pass through pelem
 896  *              (x1, y1, x2, y2 are undefined)
 897  */
 898
 899 bool
 900 PhantomElement::clipLineNormalizedCoords (double& x1, double& y1, double& x2, double& y2) const
 901 {
 902   bool accept = false;
 903
 904   switch (m_type) {
 905   case PELEM_RECTANGLE:
 906     double rect[4];
 907     rect[0] = -1.0;  rect[1] = -1.0;
 908     rect[2] = 1.0;  rect[3] = 1.0;
 909     accept = clip_rect (x1, y1, x2, y2, rect);
 910     break;
 911   case PELEM_ELLIPSE:
 912     accept = clip_circle (x1, y1, x2, y2, 0.0, 0.0, 1.0, 0.0, 0.0);
 913     break;
 914   case PELEM_TRIANGLE:
 915     accept = clip_triangle (x1, y1, x2, y2, 1.0, 1.0, true);
 916     break;
 917   case PELEM_SEGMENT:
 918     accept = clip_segment (x1, y1, x2, y2, m_u, m_v);
 919     break;
 920   case PELEM_SECTOR:
 921     accept = clip_sector (x1, y1, x2, y2, m_u, m_v);
 922     break;
 923   default:
 924     sys_error (ERR_WARNING, "Illegal pelem type %d [pelem_clip_line]", m_type);
 925     break;
 926   }
 927
 928   return(accept);
 929 }
 930
 931
 932 // METHOD IDENTIFICATION
 933 //    PhantomElement::isPointInside             Check if point is inside pelem
 934 //
 935 // SYNOPSIS
 936 //    is_point_inside (pelem, x, y, coord_type)
 937 //    double x, y               Point to see if lies in pelem
 938 //    int coord_type            Coordinate type (PELEM_COORD or PHM_COORD)
 939 //
 940 // RETURNS
 941 //    true if point lies within pelem
 942 //    false if point lies outside of pelem
 943
 944 bool
 945 PhantomElement::isPointInside (double x, double y, const CoordType coord_type) const
 946 {
 947   if (coord_type == PHM_COORD) {
 948     xform_mtx2 (m_xformPhmToObj, x, y);
 949   } else if (coord_type != PELEM_COORD) {
 950     sys_error(ERR_WARNING, "Illegal coordinate type in pelem_is_point_inside");
 951     return (false);
 952   }
 953
 954   switch (m_type) {
 955   case PELEM_RECTANGLE:
 956     if (x > 1. || x < -1. || y > 1. || y < -1.)
 957       return (false);
 958     else
 959       return (true);
 960     break;
 961   case PELEM_TRIANGLE:
 962     if (y < 0. || y > 1. - x || y > 1. + x)
 963       return (false);
 964     else
 965       return (true);
 966     break;
 967   case PELEM_ELLIPSE:
 968     if (x > 1. || x < -1. || y > 1. || y < -1.)
 969       return (false);
 970     if (x * x + y * y > 1.)             // check if inside unit circle
 971       return (false);
 972     else
 973       return (true);
 974     break;
 975
 976     // for clipping segments & sectors, must NOT scale by (1/u, 1/v)
 977     // because this destroys information about size of arc component
 978
 979   case PELEM_SEGMENT:
 980     if (x > 1. || x < -1. || y > 0.)
 981         return (false);         // clip against y > 0
 982     x *= m_u;                   // put back u & v scale
 983     y *= m_v;
 984     if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
 985       return (false);           // clip against circle, r = sqrt(@)
 986     else
 987       return (true);
 988     break;
 989   case PELEM_SECTOR:
 990       if (x > 1. || x < -1. || y > 1.)   // extent
 991       return (false);
 992       if (y > 1. - x || y > 1. + x)      // triangle
 993           return (false);                      // clip against triangle
 994       x *= m_u;                // circle: put back u & v scale
 995     y *= m_v;
 996     if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
 997         return (false);                // clip against circle
 998     else
 999       return (true);
1000   break;
1001   default:
1002     sys_error (ERR_WARNING, "Illegal pelem type in pelem_is_point_inside()");
1003     break;
1004   }
1005
1006   return (false);
1007 }
1008
1009