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