libctsim/phantom.cpp

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