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