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