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