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