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