src/NIS/NIS_Triangulated.cxx

   1 // Created on: 2007-07-17
   2 // Created by: Alexander GRIGORIEV
   3 // Copyright (c) 2007-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and / or modify it
   8 // under the terms of the GNU Lesser General Public version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <NIS_Triangulated.hxx>
  17 #include <NIS_TriangulatedDrawer.hxx>
  18 #include <NIS_InteractiveContext.hxx>
  19 #include <gp_Ax1.hxx>
  20 #include <Precision.hxx>
  21
  22 IMPLEMENT_STANDARD_HANDLE  (NIS_Triangulated, NIS_InteractiveObject)
  23 IMPLEMENT_STANDARD_RTTIEXT (NIS_Triangulated, NIS_InteractiveObject)
  24
  25 static Standard_Real aTolConf = Precision::Confusion() * 0.0001;
  26
  27 /**
  28  * Checking the given value if it is considerably greater than zero.
  29  */
  30 static inline Standard_Boolean IsPositive (const Standard_Real theVal)
  31 {
  32   return (theVal > aTolConf);
  33 }
  34
  35 /**
  36  * Compute the size in bytes of an index array.
  37  */
  38 static inline Standard_Size    NBytesInd  (const Standard_Integer nInd,
  39                                            const unsigned int     theIndType)
  40 {
  41   Standard_Size nBytes = static_cast<Standard_Size>(nInd);
  42   if (theIndType) {
  43     nBytes *= 2;
  44     if (theIndType > 1u)
  45       nBytes *= 2;
  46   }
  47   return nBytes;
  48 }
  49
  50 //=======================================================================
  51 //function : NIS_Triangulated()
  52 //purpose  : Constructor
  53 //=======================================================================
  54
  55 NIS_Triangulated::NIS_Triangulated
  56                         (const Standard_Integer                   nNodes,
  57                          const Standard_Boolean                   is2D,
  58                          const Handle(NCollection_BaseAllocator)& theAlloc)
  59   : myAlloc             (0L),
  60     myType              (Type_None),
  61     mypNodes            (0L),
  62     mypTriangles        (0L),
  63     mypLines            (0L),
  64     mypPolygons         (0L),
  65     myNNodes            (0),
  66     myNTriangles        (0),
  67     myNLineNodes        (0),
  68     myNPolygons         (0u),
  69     myIsDrawPolygons    (Standard_False),
  70     myIsCloned          (Standard_False),
  71     myIndexType         (2u),
  72     myNodeCoord         (is2D ? 2 : 3),
  73     myPolygonType       (static_cast<unsigned int>(Polygon_Default))
  74 {
  75
  76   if (theAlloc.IsNull())
  77     myAlloc = NCollection_BaseAllocator::CommonBaseAllocator().operator->();
  78   else
  79     myAlloc = theAlloc.operator->();
  80   allocateNodes (nNodes);
  81 }
  82
  83 //=======================================================================
  84 //function : Clear
  85 //purpose  : Reset all internal data members and structures
  86 //=======================================================================
  87
  88 void NIS_Triangulated::Clear ()
  89 {
  90   if (myNNodes) {
  91     myNNodes = 0;
  92     myAlloc->Free(mypNodes);
  93     mypNodes = 0L;
  94   }
  95   if (myNTriangles) {
  96     myNTriangles = 0;
  97     myAlloc->Free(mypTriangles);
  98     mypTriangles = 0L;
  99   }
 100   if (myNLineNodes) {
 101     myNLineNodes = 0;
 102     myAlloc->Free( mypLines);
 103     mypLines = 0L;
 104   }
 105   if (myNPolygons) {
 106     for (unsigned int i = 0; i < myNPolygons; i++)
 107       myAlloc->Free(mypPolygons[i]);
 108     myAlloc->Free(mypPolygons);
 109     myNPolygons = 0u;
 110     mypPolygons = 0L;
 111   }
 112   myType = Type_None;
 113   myIsDrawPolygons = Standard_False;
 114   myPolygonType = static_cast<unsigned int>(Polygon_Default);
 115   if (GetDrawer().IsNull() == Standard_False) {
 116     GetDrawer()->SetUpdated(NIS_Drawer::Draw_Normal,
 117                             NIS_Drawer::Draw_Top,
 118                             NIS_Drawer::Draw_Transparent,
 119                             NIS_Drawer::Draw_Hilighted);
 120   }
 121 }
 122
 123 //=======================================================================
 124 //function : SetPolygonsPrs
 125 //purpose  :
 126 //=======================================================================
 127
 128 void NIS_Triangulated::SetPolygonsPrs (const Standard_Integer nPolygons,
 129                                        const Standard_Integer nNodes)
 130 {
 131   if (nPolygons <= 0)
 132     myType &= ~Type_Polygons;
 133   else {
 134     myType |= Type_Polygons;
 135     if (myNPolygons) {
 136       for (unsigned int i = 0; i < myNPolygons; i++)
 137         myAlloc->Free(mypPolygons[i]);
 138       myAlloc->Free(mypPolygons);
 139     }
 140     myNPolygons = static_cast<unsigned int>(nPolygons);
 141     mypPolygons = static_cast<Standard_Integer **>
 142       (myAlloc->Allocate(sizeof(Standard_Integer *)*nPolygons));
 143     allocateNodes (nNodes);
 144   }
 145 }
 146
 147 //=======================================================================
 148 //function : SetTriangulationPrs
 149 //purpose  :
 150 //=======================================================================
 151
 152 void NIS_Triangulated::SetTriangulationPrs (const Standard_Integer nTri,
 153                                             const Standard_Integer nNodes)
 154 {
 155   if (nTri <= 0)
 156     myType &= ~Type_Triangulation;
 157   else {
 158     myType |= Type_Triangulation;
 159     if (myNTriangles)
 160       myAlloc->Free(mypTriangles);
 161     allocateNodes (nNodes);
 162
 163     myNTriangles = nTri;
 164     const Standard_Size nBytes = NBytesInd(3 * nTri, myIndexType);
 165     mypTriangles = static_cast<Standard_Integer *> (myAlloc->Allocate(nBytes));
 166   }
 167 }
 168
 169 //=======================================================================
 170 //function : SetLinePrs
 171 //purpose  :
 172 //=======================================================================
 173
 174 void NIS_Triangulated::SetLinePrs (const Standard_Integer nPoints,
 175                                    const Standard_Boolean isClosed,
 176                                    const Standard_Integer nNodes)
 177 {
 178   if (nPoints <= 0)
 179     myType &= ~(Type_Loop | Type_Line);
 180   else {
 181     myType |= Type_Line;
 182     if (isClosed)
 183       myType |= Type_Loop;
 184     if (myNLineNodes)
 185       myAlloc->Free(mypLines);
 186     myType &= ~Type_Segments;
 187     allocateNodes (nNodes);
 188
 189     myNLineNodes = nPoints;
 190     const Standard_Size nBytes = NBytesInd(nPoints, myIndexType);
 191     mypLines = static_cast<Standard_Integer *> (myAlloc->Allocate(nBytes));
 192   }
 193 }
 194
 195 //=======================================================================
 196 //function : SetSegmentPrs
 197 //purpose  :
 198 //=======================================================================
 199
 200 void NIS_Triangulated::SetSegmentPrs (const Standard_Integer nSegments,
 201                                       const Standard_Integer nNodes)
 202 {
 203   if (nSegments <= 0)
 204     myType &= ~(Type_Loop | Type_Segments);
 205   else {
 206     myType |= Type_Segments;
 207     if (myNLineNodes)
 208       myAlloc->Free(mypLines);
 209     myType &= ~(Type_Line | Type_Loop);
 210     allocateNodes (nNodes);
 211
 212     myNLineNodes = nSegments*2;
 213     const Standard_Size nBytes = NBytesInd(myNLineNodes, myIndexType);
 214     mypLines = static_cast<Standard_Integer *> (myAlloc->Allocate(nBytes));
 215   }
 216 }
 217
 218 //=======================================================================
 219 //function : ~NIS_Triangulated
 220 //purpose  : Destructor
 221 //=======================================================================
 222
 223 NIS_Triangulated::~NIS_Triangulated ()
 224 {
 225   Clear();
 226 }
 227
 228 //=======================================================================
 229 //function : DefaultDrawer
 230 //purpose  :
 231 //=======================================================================
 232
 233 NIS_Drawer * NIS_Triangulated::DefaultDrawer (NIS_Drawer * theDrawer) const
 234 {
 235   NIS_TriangulatedDrawer * aDrawer =
 236     theDrawer ? static_cast<NIS_TriangulatedDrawer *>(theDrawer)
 237               : new NIS_TriangulatedDrawer (Quantity_NOC_RED);
 238   aDrawer->myIsDrawPolygons = myIsDrawPolygons;
 239   aDrawer->myPolygonType = myPolygonType;
 240   return aDrawer;
 241 }
 242
 243 //=======================================================================
 244 //function : ComputeBox
 245 //purpose  :
 246 //=======================================================================
 247
 248 void NIS_Triangulated::ComputeBox (Bnd_B3f&                  theBox,
 249                                    const Standard_Integer    nNodes,
 250                                    const Standard_ShortReal* pNodes,
 251                                    const Standard_Integer    nCoord)
 252 {
 253   theBox.Clear();
 254   if (nNodes > 0) {
 255     Standard_ShortReal aBox[6] = {
 256       pNodes[0], pNodes[1], 0.,
 257       pNodes[0], pNodes[1], 0.
 258     };
 259     if (nCoord > 2) {
 260       aBox[2] = pNodes[2];
 261       aBox[5] = pNodes[2];
 262     }
 263     for (Standard_Integer i = 1; i < nNodes; i++) {
 264       const Standard_ShortReal * pNode = &pNodes[i * nCoord];
 265       if (aBox[0] > pNode[0])
 266         aBox[0] = pNode[0];
 267       else if (aBox[3] < pNode[0])
 268         aBox[3] = pNode[0];
 269       if (aBox[1] > pNode[1])
 270         aBox[1] = pNode[1];
 271       else if (aBox[4] < pNode[1])
 272         aBox[4] = pNode[1];
 273       if (nCoord > 2) {
 274         if (aBox[2] > pNode[2])
 275           aBox[2] = pNode[2];
 276         else if (aBox[5] < pNode[2])
 277           aBox[5] = pNode[2];
 278       }
 279     }
 280     theBox.Add (gp_XYZ (Standard_Real(aBox[0]),
 281                         Standard_Real(aBox[1]),
 282                         Standard_Real(aBox[2])));
 283     theBox.Add (gp_XYZ (Standard_Real(aBox[3]),
 284                         Standard_Real(aBox[4]),
 285                         Standard_Real(aBox[5])));
 286   }
 287 }
 288
 289 //=======================================================================
 290 //function : computeBox
 291 //purpose  :
 292 //=======================================================================
 293
 294 void NIS_Triangulated::computeBox ()
 295 {
 296   ComputeBox (myBox, myNNodes, mypNodes, myNodeCoord);
 297 }
 298
 299 //=======================================================================
 300 //function : SetNode
 301 //purpose  :
 302 //=======================================================================
 303
 304 void NIS_Triangulated::SetNode       (const Standard_Integer  ind,
 305                                       const gp_XYZ&           thePnt)
 306 {
 307   if (ind >= myNNodes)
 308     Standard_OutOfRange::Raise ("NIS_Triangulated::SetNode");
 309   Standard_ShortReal * pNode = &mypNodes[myNodeCoord * ind];
 310   pNode[0] = Standard_ShortReal(thePnt.X());
 311   pNode[1] = Standard_ShortReal(thePnt.Y());
 312   if (myNodeCoord > 2)
 313     pNode[2] = Standard_ShortReal(thePnt.Z());
 314   setIsUpdateBox(Standard_True);
 315 }
 316
 317 //=======================================================================
 318 //function : SetNode
 319 //purpose  :
 320 //=======================================================================
 321
 322 void NIS_Triangulated::SetNode       (const Standard_Integer  ind,
 323                                       const gp_XY&            thePnt)
 324 {
 325   if (ind >= myNNodes)
 326     Standard_OutOfRange::Raise ("NIS_Triangulated::SetNode");
 327   Standard_ShortReal * pNode = &mypNodes[myNodeCoord * ind];
 328   pNode[0] = Standard_ShortReal(thePnt.X());
 329   pNode[1] = Standard_ShortReal(thePnt.Y());
 330   if (myNodeCoord > 2)
 331     pNode[2] = 0.f;
 332   setIsUpdateBox(Standard_True);
 333 }
 334
 335 //=======================================================================
 336 //function : SetTriangle
 337 //purpose  :
 338 //=======================================================================
 339
 340 void NIS_Triangulated::SetTriangle   (const Standard_Integer  ind,
 341                                       const Standard_Integer  iNode0,
 342                                       const Standard_Integer  iNode1,
 343                                       const Standard_Integer  iNode2)
 344 {
 345   if (ind >= myNTriangles)
 346     Standard_OutOfRange::Raise ("NIS_Triangulated::SetTriangle");
 347   switch (myIndexType) {
 348     case 0: // 8bit
 349     {
 350       unsigned char * pTri =
 351         reinterpret_cast<unsigned char *>(mypTriangles) + (3 * ind);
 352       pTri[0] = static_cast<unsigned char>(iNode0);
 353       pTri[1] = static_cast<unsigned char>(iNode1);
 354       pTri[2] = static_cast<unsigned char>(iNode2);
 355     }
 356     break;
 357     case 1: // 16bit
 358     {
 359       unsigned short * pTri =
 360         reinterpret_cast<unsigned short *>(mypTriangles) + (3 * ind);
 361       pTri[0] = static_cast<unsigned short>(iNode0);
 362       pTri[1] = static_cast<unsigned short>(iNode1);
 363       pTri[2] = static_cast<unsigned short>(iNode2);
 364     }
 365     break;
 366     default: // 32bit
 367     {
 368       Standard_Integer * pTri = &mypTriangles[3*ind];
 369       pTri[0] = iNode0;
 370       pTri[1] = iNode1;
 371       pTri[2] = iNode2;
 372     }
 373   }
 374 }
 375
 376 //=======================================================================
 377 //function : SetLineNode
 378 //purpose  :
 379 //=======================================================================
 380
 381 void NIS_Triangulated::SetLineNode      (const Standard_Integer ind,
 382                                          const Standard_Integer iNode)
 383 {
 384   if (ind >= myNLineNodes)
 385     Standard_OutOfRange::Raise ("NIS_Triangulated::SetTriangle");
 386   switch (myIndexType) {
 387     case 0: // 8bit
 388     {
 389       unsigned char * pInd =
 390         reinterpret_cast<unsigned char *>(mypLines) + ind;
 391       pInd[0] = static_cast<unsigned char>(iNode);
 392     }
 393     break;
 394     case 1: // 16bit
 395     {
 396       unsigned short * pInd =
 397         reinterpret_cast<unsigned short *>(mypLines) + ind;
 398       pInd[0] = static_cast<unsigned short>(iNode);
 399     }
 400     break;
 401     default: // 32bit
 402       mypLines[ind] = iNode;
 403   }
 404 }
 405
 406 //=======================================================================
 407 //function : SetPolygonNode
 408 //purpose  :
 409 //=======================================================================
 410
 411 void NIS_Triangulated::SetPolygonNode         (const Standard_Integer indPoly,
 412                                                const Standard_Integer ind,
 413                                                const Standard_Integer iNode)
 414 {
 415   if (indPoly >= static_cast<Standard_Integer>(myNPolygons))
 416     Standard_OutOfRange::Raise ("NIS_Triangulated::SetPolygonNode");
 417
 418   Standard_Integer * aPoly  = mypPolygons[indPoly];
 419   switch (myIndexType) {
 420   case 0: // 8bit
 421   {
 422     unsigned char aNNode =  * (reinterpret_cast<unsigned char *>(aPoly));
 423     if (static_cast<unsigned char>(ind) >= aNNode)
 424       Standard_OutOfRange::Raise ("NIS_Triangulated::SetPolygonNode");
 425
 426     unsigned char * pInd =
 427       reinterpret_cast<unsigned char *>(aPoly) + ind + 1;
 428     pInd[0] = static_cast<unsigned char>(iNode);
 429   }
 430   break;
 431   case 1: // 16bit
 432   {
 433     unsigned short aNNode =  * (reinterpret_cast<unsigned short *>(aPoly));
 434     if (static_cast<unsigned short>(ind) >= aNNode)
 435       Standard_OutOfRange::Raise ("NIS_Triangulated::SetPolygonNode");
 436
 437     unsigned short * pInd =
 438       reinterpret_cast<unsigned short *>(aPoly) + ind + 1;
 439     pInd[0] = static_cast<unsigned short>(iNode);
 440   }
 441   break;
 442   default: // 32bit
 443     if (ind >= aPoly[0])
 444       Standard_OutOfRange::Raise ("NIS_Triangulated::SetPolygonNode");
 445     aPoly[ind + 1] = iNode;
 446   }
 447 }
 448
 449 //=======================================================================
 450 //function : NPolygonNodes
 451 //purpose  :
 452 //=======================================================================
 453 Standard_Integer NIS_Triangulated::NPolygonNodes
 454                                          (const Standard_Integer indPoly)const
 455 {
 456   Standard_Integer aResult(0);
 457   if (indPoly >= static_cast<Standard_Integer>(myNPolygons))
 458     Standard_OutOfRange::Raise ("NIS_Triangulated::PolygonNode");
 459   Standard_Integer * aPoly  = mypPolygons[indPoly];
 460   switch (myIndexType) {
 461     case 0: // 8bit
 462     {
 463       unsigned char aNNode =  * (reinterpret_cast<unsigned char *>(aPoly));
 464       aResult = static_cast<Standard_Integer>(aNNode);
 465     }
 466     break;
 467     case 1: // 16bit
 468     {
 469       unsigned short aNNode =  * (reinterpret_cast<unsigned short *>(aPoly));
 470       aResult = static_cast<Standard_Integer>(aNNode);
 471     }
 472     break;
 473     default: // 32bit
 474     {
 475       aResult = aPoly[0];
 476     }
 477   }
 478   return aResult;
 479 }
 480
 481 //=======================================================================
 482 //function : PolygonNode
 483 //purpose  :
 484 //=======================================================================
 485
 486 Standard_Integer NIS_Triangulated::PolygonNode
 487                                             (const Standard_Integer indPoly,
 488                                              const Standard_Integer ind)const
 489 {
 490   Standard_Integer aResult(-1);
 491   if (indPoly >= static_cast<Standard_Integer>(myNPolygons))
 492     Standard_OutOfRange::Raise ("NIS_Triangulated::PolygonNode");
 493   const Standard_Integer * aPoly  = mypPolygons[indPoly];
 494   switch (myIndexType) {
 495     case 0: // 8bit
 496     {
 497       const unsigned char * pInd =
 498         reinterpret_cast<const unsigned char *>(aPoly);
 499       if (static_cast<unsigned char>(ind) >= pInd[0])
 500         Standard_OutOfRange::Raise ("NIS_Triangulated::PolygonNode");
 501
 502       aResult = static_cast<Standard_Integer>(pInd[ind + 1]);
 503     }
 504     break;
 505     case 1: // 16bit
 506     {
 507       const unsigned short * pInd =
 508         reinterpret_cast<const unsigned short *>(aPoly);
 509       if (static_cast<unsigned short>(ind) >= pInd[0])
 510         Standard_OutOfRange::Raise ("NIS_Triangulated::PolygonNode");
 511
 512       aResult = static_cast<Standard_Integer>(pInd[ind + 1]);
 513     }
 514     break;
 515     default: // 32bit
 516     {
 517       if (ind >= aPoly[0])
 518         Standard_OutOfRange::Raise ("NIS_Triangulated::PolygonNode");
 519       aResult = aPoly[ind + 1];
 520     }
 521   }
 522   return aResult;
 523 }
 524
 525 //=======================================================================
 526 //function : SetPolygon
 527 //purpose  :
 528 //=======================================================================
 529
 530 void NIS_Triangulated::SetPolygon (const Standard_Integer  ind,
 531                                    const Standard_Integer  theSz)
 532 {
 533   if (ind >= static_cast<Standard_Integer>(myNPolygons))
 534     Standard_OutOfRange::Raise ("NIS_Triangulated::SetPolygon");
 535   switch (myIndexType) {
 536     case 0: // 8bit
 537     {
 538       unsigned char * anArray  = static_cast <unsigned char *>
 539         (myAlloc->Allocate (sizeof(unsigned char) * (theSz+1)));
 540       mypPolygons[ind] = reinterpret_cast<Standard_Integer *> (anArray);
 541       anArray[0] = static_cast<unsigned char>(theSz);
 542     }
 543     break;
 544     case 1: // 16bit
 545     {
 546       unsigned short * anArray  = static_cast <unsigned short *>
 547         (myAlloc->Allocate (sizeof(unsigned short) * (theSz+1)));
 548       mypPolygons[ind] = reinterpret_cast<Standard_Integer *> (anArray);
 549       anArray[0] = static_cast<unsigned short>(theSz);
 550     }
 551     break;
 552     default: // 32bit
 553     {
 554       Standard_Integer * anArray  = static_cast <Standard_Integer *>
 555         (myAlloc->Allocate (sizeof(Standard_Integer) * (theSz+1)));
 556       mypPolygons[ind] = anArray;
 557       anArray[0] = theSz;
 558     }
 559   }
 560 }
 561
 562 //=======================================================================
 563 //function : SetDrawPolygons
 564 //purpose  :
 565 //=======================================================================
 566
 567 void NIS_Triangulated::SetDrawPolygons (const Standard_Boolean isDrawPolygons)
 568 {
 569   if (GetDrawer().IsNull())
 570     myIsDrawPolygons = isDrawPolygons;
 571   else {
 572     if (myIsDrawPolygons != isDrawPolygons) {
 573       const Handle(NIS_TriangulatedDrawer) aDrawer =
 574         static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 575       aDrawer->Assign (GetDrawer());
 576       aDrawer->myIsDrawPolygons = isDrawPolygons;
 577       SetDrawer (aDrawer);
 578       myIsDrawPolygons = isDrawPolygons;
 579     }
 580   }
 581 }
 582
 583 //=======================================================================
 584 //function : SetPolygonType
 585 //purpose  : Set the type of polygon rendering
 586 //=======================================================================
 587
 588 void NIS_Triangulated::SetPolygonType
 589                                 (const NIS_Triangulated::PolygonType theType)
 590 {
 591   if (GetDrawer().IsNull())
 592     myPolygonType = static_cast<unsigned int>(theType);
 593   else {
 594     if (myPolygonType != static_cast<unsigned int>(theType)) {
 595       const Handle(NIS_TriangulatedDrawer) aDrawer =
 596         static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 597       aDrawer->Assign (GetDrawer());
 598       aDrawer->myPolygonType = theType;
 599       SetDrawer (aDrawer);
 600       myPolygonType = static_cast<unsigned int>(theType);
 601     }
 602   }
 603 }
 604
 605 //=======================================================================
 606 //function : SetColor
 607 //purpose  : Set the normal color for presentation.
 608 //=======================================================================
 609
 610 void NIS_Triangulated::SetColor (const Quantity_Color&  theColor)
 611 {
 612   const Handle(NIS_TriangulatedDrawer) aDrawer =
 613     static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 614   aDrawer->Assign (GetDrawer());
 615   aDrawer->myColor[NIS_Drawer::Draw_Normal] = theColor;
 616   aDrawer->myColor[NIS_Drawer::Draw_Top] = theColor;
 617   aDrawer->myColor[NIS_Drawer::Draw_Transparent] = theColor;
 618   SetDrawer (aDrawer);
 619 }
 620
 621 //=======================================================================
 622 //function : GetColor
 623 //purpose  : Get Normal, Transparent or Hilighted color of the presentation.
 624 //=======================================================================
 625
 626 Quantity_Color NIS_Triangulated::GetColor
 627                         (const NIS_Drawer::DrawType theDrawType) const
 628 {
 629   Handle(NIS_TriangulatedDrawer) aDrawer =
 630     Handle(NIS_TriangulatedDrawer)::DownCast( GetDrawer() );
 631   if (aDrawer.IsNull() == Standard_False)
 632   {
 633     return aDrawer->myColor[theDrawType];
 634   }
 635   return Quantity_Color(); // return null color object
 636 }
 637
 638 //=======================================================================
 639 //function : SetHilightColor
 640 //purpose  : Set the color for hilighted presentation.
 641 //=======================================================================
 642
 643 void NIS_Triangulated::SetHilightColor (const Quantity_Color&  theColor)
 644 {
 645   const Handle(NIS_TriangulatedDrawer) aDrawer =
 646     static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 647   aDrawer->Assign (GetDrawer());
 648   aDrawer->myColor[NIS_Drawer::Draw_Hilighted] = theColor;
 649   SetDrawer (aDrawer);
 650 }
 651
 652 //=======================================================================
 653 //function : SetDynHilightColor
 654 //purpose  : Set the color for dynamic hilight presentation.
 655 //=======================================================================
 656
 657 void NIS_Triangulated::SetDynHilightColor(const Quantity_Color&  theColor)
 658 {
 659   const Handle(NIS_TriangulatedDrawer) aDrawer =
 660     static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 661   aDrawer->Assign (GetDrawer());
 662   aDrawer->myColor[NIS_Drawer::Draw_DynHilighted] = theColor;
 663   SetDrawer (aDrawer);
 664 }
 665
 666 //=======================================================================
 667 //function : SetLineWidth
 668 //purpose  : Set the width of line presentations in pixels.
 669 //=======================================================================
 670
 671 void NIS_Triangulated::SetLineWidth (const Standard_Real    theWidth)
 672 {
 673   const Handle(NIS_TriangulatedDrawer) aDrawer =
 674     static_cast<NIS_TriangulatedDrawer *>(DefaultDrawer(0L));
 675   aDrawer->Assign (GetDrawer());
 676   aDrawer->myLineWidth = (Standard_ShortReal) theWidth;
 677   SetDrawer (aDrawer);
 678 }
 679
 680 //=======================================================================
 681 //function : Clone
 682 //purpose  :
 683 //=======================================================================
 684
 685 void NIS_Triangulated::Clone (const Handle_NCollection_BaseAllocator& theAlloc,
 686                               Handle_NIS_InteractiveObject& theDest) const
 687 {
 688   Handle(NIS_Triangulated) aNewObj;
 689   if (theDest.IsNull()) {
 690     aNewObj = new (theAlloc) NIS_Triangulated(myNNodes, myNodeCoord == 2,
 691                                               theAlloc);
 692     aNewObj->myIsCloned = Standard_True;
 693     theDest = aNewObj;
 694   } else {
 695     aNewObj = reinterpret_cast<NIS_Triangulated*> (theDest.operator->());
 696     aNewObj->myAlloc = theAlloc.operator->();
 697     aNewObj->myNNodes = 0;
 698     aNewObj->allocateNodes(myNNodes);
 699   }
 700   NIS_InteractiveObject::Clone(theAlloc, theDest);
 701   if (myNNodes > 0)
 702   {
 703     // copy nodes
 704     memcpy(aNewObj->mypNodes, mypNodes,
 705            myNNodes * myNodeCoord * sizeof(Standard_ShortReal));
 706     // copy triangles
 707     aNewObj->myNTriangles = myNTriangles;
 708     if (myNTriangles) {
 709       const Standard_Size nBytes = NBytesInd(3 * myNTriangles, myIndexType);
 710       aNewObj->mypTriangles = static_cast<Standard_Integer *>
 711         (theAlloc->Allocate(nBytes));
 712       memcpy(aNewObj->mypTriangles, mypTriangles, nBytes);
 713     }
 714     // copy lines/segments
 715     aNewObj->myNLineNodes = myNLineNodes;
 716     if (myNLineNodes) {
 717       const Standard_Size nBytes = NBytesInd(myNLineNodes, myIndexType);
 718       aNewObj->mypLines = static_cast<Standard_Integer *>
 719         (theAlloc->Allocate(nBytes));
 720       memcpy(aNewObj->mypLines, mypLines, nBytes);
 721     }
 722     // copy polygons
 723     aNewObj->myNPolygons = myNPolygons;
 724     if (myNPolygons) {
 725       const Standard_Size nBytes = sizeof(Standard_Integer *)*myNPolygons;
 726       aNewObj->mypPolygons = static_cast<Standard_Integer **>
 727         (theAlloc->Allocate(nBytes));
 728       for (unsigned int i = 0; i < myNPolygons; i++) {
 729         const Standard_Integer nNodes = NPolygonNodes(i);
 730         const Standard_Size nBytes = NBytesInd(nNodes+1, myIndexType);
 731         aNewObj->mypPolygons[i] = static_cast <Standard_Integer *>
 732           (theAlloc->Allocate (nBytes));
 733         memcpy(aNewObj->mypPolygons[i], mypPolygons[i], nBytes);
 734       }
 735     }
 736   }
 737   aNewObj->myType = myType;
 738   aNewObj->myIsDrawPolygons = myIsDrawPolygons;
 739   aNewObj->myIndexType = myIndexType;
 740   aNewObj->myPolygonType = myPolygonType;
 741 }
 742
 743 //=======================================================================
 744 //function : Delete
 745 //purpose  :
 746 //=======================================================================
 747
 748 void NIS_Triangulated::Delete () const
 749 {
 750   if (myIsCloned == Standard_False)
 751     Standard_Transient::Delete();
 752   else {
 753     // Call the destructor and then release the memory occupied by the instance.
 754     // This is required when the NIS_Triangulated instance is allocated in
 755     // the same allocator as its internal arrays.
 756     NIS_Triangulated* pThis = const_cast<NIS_Triangulated*>(this);
 757     pThis->~NIS_Triangulated();
 758     myAlloc->Free(pThis);
 759   }
 760 }
 761
 762 //=======================================================================
 763 //function : Intersect
 764 //purpose  :
 765 //=======================================================================
 766
 767 Standard_Boolean NIS_Triangulated::Intersect
 768                                         (const Bnd_B3f&         theBox,
 769                                          const gp_Trsf&         theTrf,
 770                                          const Standard_Boolean isFullIn) const
 771 {
 772   Standard_Boolean aResult (isFullIn);
 773
 774   if ((myType & Type_Triangulation) && myIsDrawPolygons == Standard_False) {
 775     unsigned int iNode = 0;
 776     for (; iNode < myNNodes * myNodeCoord; iNode += myNodeCoord)
 777     {
 778       gp_XYZ aPnt (static_cast<Standard_Real>(mypNodes[iNode+0]),
 779                    static_cast<Standard_Real>(mypNodes[iNode+1]), 0.);
 780       if (myNodeCoord > 2)
 781         aPnt.SetZ (static_cast<Standard_Real>(mypNodes[iNode+2]));
 782       theTrf.Transforms (aPnt);
 783       if (theBox.IsOut (aPnt)) {
 784         if (isFullIn) {
 785           aResult = Standard_False;
 786           break;
 787         }
 788       } else {
 789         if (isFullIn == Standard_False) {
 790           aResult = Standard_True;
 791           break;
 792         }
 793       }
 794     }
 795   }
 796   if (aResult == isFullIn) {
 797     if (myType & Type_Segments) {
 798       for (Standard_Integer i = 0; i < myNLineNodes; i+=2) {
 799         const gp_Pnt aPnt[2] = {
 800           nodeAtInd(mypLines, i+0).Transformed(theTrf),
 801           nodeAtInd(mypLines, i+1).Transformed(theTrf)
 802         };
 803         if (isFullIn) {
 804           if (seg_box_included (theBox, aPnt) == 0) {
 805             aResult = Standard_False;
 806             break;
 807           }
 808         } else {
 809           if (seg_box_intersect (theBox, aPnt)) {
 810             aResult = Standard_True;
 811             break;
 812           }
 813         }
 814       }
 815     } else if (myType & Type_Line) {
 816       for (Standard_Integer i = 1; i < myNLineNodes; i++) {
 817         const gp_Pnt aPnt[2] = {
 818           nodeAtInd(mypLines, i-1).Transformed(theTrf),
 819           nodeAtInd(mypLines, i+0).Transformed(theTrf)
 820         };
 821         if (isFullIn) {
 822           if (seg_box_included (theBox, aPnt) == 0) {
 823             aResult = Standard_False;
 824             break;
 825           }
 826         } else {
 827           if (seg_box_intersect (theBox, aPnt)) {
 828             aResult = Standard_True;
 829             break;
 830           }
 831         }
 832       }
 833       if (aResult == isFullIn && (myType & Type_Loop)) {
 834         const gp_Pnt aPntLast[2] = {
 835           nodeAtInd(mypLines, myNLineNodes-1).Transformed(theTrf),
 836           nodeAtInd(mypLines, 0).Transformed(theTrf)
 837         };
 838         if (isFullIn) {
 839           if (seg_box_included (theBox, aPntLast) == 0)
 840             aResult = Standard_False;
 841         } else {
 842           if (seg_box_intersect (theBox, aPntLast))
 843             aResult = Standard_True;
 844         }
 845       }
 846     } else if ((myType & Type_Polygons) && myIsDrawPolygons) {
 847       for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
 848         const Standard_Integer * aPoly = mypPolygons[iPoly];
 849         const Standard_Integer nNodes = NPolygonNodes(iPoly);
 850         for (Standard_Integer i = 1; i < nNodes; i++) {
 851           // index is incremented by 1 for the head number in the array
 852           const gp_Pnt aPnt[2] = {
 853             nodeAtInd(aPoly, i+0).Transformed(theTrf),
 854             nodeAtInd(aPoly, i+1).Transformed(theTrf)
 855           };
 856           if (isFullIn) {
 857             if (seg_box_included (theBox, aPnt) == 0) {
 858               aResult = Standard_False;
 859               break;
 860             }
 861           } else {
 862             if (seg_box_intersect (theBox, aPnt)) {
 863               aResult = Standard_True;
 864               break;
 865             }
 866           }
 867         }
 868         if (aResult == isFullIn) {
 869           const gp_Pnt aPntLast[2] = {
 870             nodeAtInd(aPoly, nNodes).Transformed(theTrf),
 871             nodeAtInd(aPoly, 1).Transformed(theTrf)
 872           };
 873           if (isFullIn) {
 874             if (seg_box_included (theBox, aPntLast) == 0)
 875               aResult = Standard_False;
 876           } else {
 877             if (seg_box_intersect (theBox, aPntLast))
 878               aResult = Standard_True;
 879           }
 880         }
 881       }
 882     }
 883   }
 884   return aResult;
 885 }
 886
 887 //=======================================================================
 888 //function : Intersect
 889 //purpose  :
 890 //=======================================================================
 891
 892 Standard_Real NIS_Triangulated::Intersect (const gp_Ax1&       theAxis,
 893                                            const Standard_Real theOver) const
 894 {
 895   Standard_Real aResult (RealLast());
 896   Standard_Real start[3], dir[3];
 897   theAxis.Location().Coord(start[0], start[1], start[2]);
 898   theAxis.Direction().Coord(dir[0], dir[1], dir[2]);
 899   double anInter;
 900
 901   if ((myType & Type_Triangulation) && (myIsDrawPolygons == Standard_False))
 902     for (Standard_Integer i = 0; i < myNTriangles; i++) {
 903       Standard_Boolean isIntersect(Standard_False);
 904       if (myIndexType == 0) {
 905         const unsigned char * pTri =
 906           reinterpret_cast<unsigned char *>(mypTriangles) + (3 * i);
 907         if (myNodeCoord > 2)
 908           isIntersect = tri_line_intersect (start, dir,
 909                                             &mypNodes[myNodeCoord * pTri[0]],
 910                                             &mypNodes[myNodeCoord * pTri[1]],
 911                                             &mypNodes[myNodeCoord * pTri[2]],
 912                                             &anInter);
 913         else
 914           isIntersect = tri2d_line_intersect (start, dir,
 915                                               &mypNodes[myNodeCoord * pTri[0]],
 916                                               &mypNodes[myNodeCoord * pTri[1]],
 917                                               &mypNodes[myNodeCoord * pTri[2]],
 918                                               &anInter);
 919       } else if (myIndexType == 1) {
 920         const unsigned short * pTri =
 921           reinterpret_cast<unsigned short *>(mypTriangles) + (3 * i);
 922         if (myNodeCoord > 2)
 923           isIntersect = tri_line_intersect (start, dir,
 924                                             &mypNodes[myNodeCoord * pTri[0]],
 925                                             &mypNodes[myNodeCoord * pTri[1]],
 926                                             &mypNodes[myNodeCoord * pTri[2]],
 927                                             &anInter);
 928         else
 929           isIntersect = tri2d_line_intersect (start, dir,
 930                                               &mypNodes[myNodeCoord * pTri[0]],
 931                                               &mypNodes[myNodeCoord * pTri[1]],
 932                                               &mypNodes[myNodeCoord * pTri[2]],
 933                                               &anInter);
 934       } else {
 935         const Standard_Integer * pTri = &mypTriangles[3 * i];
 936         if (myNodeCoord > 2)
 937           isIntersect = tri_line_intersect (start, dir,
 938                                             &mypNodes[myNodeCoord * pTri[0]],
 939                                             &mypNodes[myNodeCoord * pTri[1]],
 940                                             &mypNodes[myNodeCoord * pTri[2]],
 941                                             &anInter);
 942         else
 943           isIntersect = tri2d_line_intersect (start, dir,
 944                                               &mypNodes[myNodeCoord * pTri[0]],
 945                                               &mypNodes[myNodeCoord * pTri[1]],
 946                                               &mypNodes[myNodeCoord * pTri[2]],
 947                                               &anInter);
 948       }
 949       if (isIntersect && anInter < aResult)
 950         aResult = anInter;
 951     }
 952
 953   const Standard_Real anOver2 = theOver*theOver;
 954   if (myType & Type_Segments) {
 955     Standard_Integer i = 0;
 956     if (myNodeCoord > 2)
 957       for (; i < myNLineNodes; i+=2) {
 958         if (seg_line_intersect (theAxis.Location().XYZ(),
 959                                 theAxis.Direction().XYZ(), anOver2,
 960                                 nodeArrAtInd(mypLines, i+0),
 961                                 nodeArrAtInd(mypLines, i+1),
 962                                 &anInter))
 963           if (anInter < aResult)
 964             aResult = anInter;
 965       }
 966     else
 967       for (; i < myNLineNodes; i+=2) {
 968         if (seg2d_line_intersect (theAxis.Location().XYZ(),
 969                                   theAxis.Direction().XYZ(), anOver2,
 970                                   nodeArrAtInd(mypLines, i+0),
 971                                   nodeArrAtInd(mypLines, i+1),
 972                                   &anInter))
 973           if (anInter < aResult)
 974             aResult = anInter;
 975       }
 976   } else if (myType & Type_Line) {
 977     Standard_Integer i = 1;
 978     if (myNodeCoord > 2) {
 979       for (; i < myNLineNodes; i++) {
 980         if (seg_line_intersect (theAxis.Location().XYZ(),
 981                                 theAxis.Direction().XYZ(), anOver2,
 982                                 nodeArrAtInd(mypLines, i-1),
 983                                 nodeArrAtInd(mypLines, i-0),
 984                                 &anInter))
 985           if (anInter < aResult)
 986             aResult = anInter;
 987       }
 988       if (myType & Type_Loop)
 989         if (seg_line_intersect (theAxis.Location().XYZ(),
 990                                 theAxis.Direction().XYZ(), anOver2,
 991                                 nodeArrAtInd(mypLines, myNLineNodes-1),
 992                                 nodeArrAtInd(mypLines, 0),
 993                                 &anInter))
 994           if (anInter < aResult)
 995             aResult = anInter;
 996     } else {
 997       for (; i < myNLineNodes; i++) {
 998         if (seg2d_line_intersect (theAxis.Location().XYZ(),
 999                                   theAxis.Direction().XYZ(), anOver2,
1000                                   nodeArrAtInd(mypLines, i-1),
1001                                   nodeArrAtInd(mypLines, i-0),
1002                                   &anInter))
1003           if (anInter < aResult)
1004             aResult = anInter;
1005       }
1006       if (myType & Type_Loop)
1007         if (seg2d_line_intersect (theAxis.Location().XYZ(),
1008                                   theAxis.Direction().XYZ(), anOver2,
1009                                   nodeArrAtInd(mypLines, myNLineNodes-1),
1010                                   nodeArrAtInd(mypLines, 0),
1011                                   &anInter))
1012           if (anInter < aResult)
1013             aResult = anInter;
1014     }
1015   }
1016
1017   if ((myType & Type_Polygons) && myIsDrawPolygons) {
1018     for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
1019       const Standard_Integer * aPoly = mypPolygons[iPoly];
1020       const Standard_Integer nNodes = NPolygonNodes(iPoly);
1021       Standard_Integer i = 1;
1022       if (myNodeCoord > 2) {
1023         for (; i < nNodes; i++) {
1024           // Node index is incremented for the head of polygon indice array
1025           if (seg_line_intersect (theAxis.Location().XYZ(),
1026                                   theAxis.Direction().XYZ(), anOver2,
1027                                   nodeArrAtInd(aPoly, i+0),
1028                                   nodeArrAtInd(aPoly, i+1),
1029                                   &anInter))
1030             if (anInter < aResult)
1031               aResult = anInter;
1032         }
1033         if (seg_line_intersect (theAxis.Location().XYZ(),
1034                                 theAxis.Direction().XYZ(), anOver2,
1035                                 nodeArrAtInd(aPoly, nNodes),
1036                                 nodeArrAtInd(aPoly, 1),
1037                                 &anInter))
1038           if (anInter < aResult)
1039             aResult = anInter;
1040       } else {
1041         for (; i < nNodes; i++) {
1042           // Node index is incremented for the head of polygon indice array
1043           if (seg2d_line_intersect (theAxis.Location().XYZ(),
1044                                     theAxis.Direction().XYZ(), anOver2,
1045                                     nodeArrAtInd(aPoly, i+0),
1046                                     nodeArrAtInd(aPoly, i+1),
1047                                     &anInter))
1048             if (anInter < aResult)
1049               aResult = anInter;
1050         }
1051         if (seg2d_line_intersect (theAxis.Location().XYZ(),
1052                                   theAxis.Direction().XYZ(), anOver2,
1053                                   nodeArrAtInd(aPoly, nNodes),
1054                                   nodeArrAtInd(aPoly, 1),
1055                                   &anInter))
1056           if (anInter < aResult)
1057             aResult = anInter;
1058       }
1059     }
1060   }
1061   return aResult;
1062 }
1063
1064 //=======================================================================
1065 //function : Intersect
1066 //purpose  :
1067 //=======================================================================
1068 Standard_Boolean NIS_Triangulated::Intersect
1069                     (const NCollection_List<gp_XY> &thePolygon,
1070                      const gp_Trsf                 &theTrf,
1071                      const Standard_Boolean         isFullIn) const
1072 {
1073   Standard_Boolean aResult (isFullIn);
1074
1075   if ((myType & Type_Triangulation) && myIsDrawPolygons == Standard_False) {
1076     unsigned int iNode = 0;
1077     for (; iNode < myNNodes * myNodeCoord; iNode += myNodeCoord)
1078     {
1079       gp_XYZ aPnt (static_cast<Standard_Real>(mypNodes[iNode+0]),
1080                    static_cast<Standard_Real>(mypNodes[iNode+1]), 0.);
1081       if (myNodeCoord > 2)
1082         aPnt.SetZ (static_cast<Standard_Real>(mypNodes[iNode+2]));
1083       theTrf.Transforms (aPnt);
1084
1085       gp_XY aP2d(aPnt.X(), aPnt.Y());
1086
1087       if (!NIS_Triangulated::IsIn(thePolygon, aP2d)) {
1088         if (isFullIn) {
1089           aResult = Standard_False;
1090           break;
1091         }
1092       } else {
1093         if (isFullIn == Standard_False) {
1094           aResult = Standard_True;
1095           break;
1096         }
1097       }
1098     }
1099   }
1100   if (aResult == isFullIn) {
1101     if (myType & Type_Segments) {
1102       for (Standard_Integer i = 0; i < myNLineNodes; i+=2) {
1103         const gp_Pnt aPnt[2] = {
1104           nodeAtInd(mypLines, i+0).Transformed(theTrf),
1105           nodeAtInd(mypLines, i+1).Transformed(theTrf)
1106         };
1107         const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1108                                 gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1109
1110         if (isFullIn) {
1111           if (seg_polygon_included (thePolygon, aP2d) == 0) {
1112             aResult = Standard_False;
1113             break;
1114           }
1115         } else {
1116           if (seg_polygon_intersect (thePolygon, aP2d)) {
1117             aResult = Standard_True;
1118             break;
1119           }
1120         }
1121       }
1122     } else if (myType & Type_Line) {
1123       for (Standard_Integer i = 1; i < myNLineNodes; i++) {
1124         const gp_Pnt aPnt[2] = {
1125           nodeAtInd(mypLines, i-1).Transformed(theTrf),
1126           nodeAtInd(mypLines, i+0).Transformed(theTrf)
1127         };
1128         const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1129                                 gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1130         if (isFullIn) {
1131           if (seg_polygon_included (thePolygon, aP2d) == 0) {
1132             aResult = Standard_False;
1133             break;
1134           }
1135         } else {
1136           if (seg_polygon_intersect (thePolygon, aP2d)) {
1137             aResult = Standard_True;
1138             break;
1139           }
1140         }
1141       }
1142       if (aResult == isFullIn && (myType & Type_Loop)) {
1143         const gp_Pnt aPntLast[2] = {
1144           nodeAtInd(mypLines, myNLineNodes-1).Transformed(theTrf),
1145           nodeAtInd(mypLines, 0).Transformed(theTrf)
1146         };
1147         const gp_XY aP2dLast[2] = { gp_XY(aPntLast[0].X(), aPntLast[0].Y()),
1148                                     gp_XY(aPntLast[1].X(), aPntLast[1].Y()) };
1149
1150         if (isFullIn) {
1151           if (seg_polygon_included (thePolygon, aP2dLast) == 0)
1152             aResult = Standard_False;
1153         } else {
1154           if (seg_polygon_intersect (thePolygon, aP2dLast))
1155             aResult = Standard_True;
1156         }
1157       }
1158     } else if ((myType & Type_Polygons) && myIsDrawPolygons) {
1159       for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
1160         const Standard_Integer * aPoly = mypPolygons[iPoly];
1161         const Standard_Integer nNodes = NPolygonNodes(iPoly);
1162         for (Standard_Integer i = 1; i < nNodes; i++) {
1163           const gp_Pnt aPnt[2] = {
1164             nodeAtInd(aPoly, i+0).Transformed(theTrf),
1165             nodeAtInd(aPoly, i+1).Transformed(theTrf)
1166           };
1167           const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1168                                   gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1169           if (isFullIn) {
1170             if (seg_polygon_included (thePolygon, aP2d) == 0) {
1171               aResult = Standard_False;
1172               break;
1173             }
1174           } else {
1175             if (seg_polygon_intersect (thePolygon, aP2d)) {
1176               aResult = Standard_True;
1177               break;
1178             }
1179           }
1180         }
1181         if (aResult == isFullIn) {
1182           const gp_Pnt aPntLast[2] = {
1183             nodeAtInd(aPoly, nNodes).Transformed(theTrf),
1184             nodeAtInd(aPoly, 1).Transformed(theTrf)
1185           };
1186           const gp_XY aP2dLast[2] = { gp_XY(aPntLast[0].X(), aPntLast[0].Y()),
1187                                       gp_XY(aPntLast[1].X(), aPntLast[1].Y()) };
1188           if (isFullIn) {
1189             if (seg_polygon_included (thePolygon, aP2dLast) == 0)
1190               aResult = Standard_False;
1191           } else {
1192             if (seg_polygon_intersect (thePolygon, aP2dLast))
1193               aResult = Standard_True;
1194           }
1195         }
1196       }
1197     }
1198   }
1199   return aResult;
1200 }
1201
1202 /* =====================================================================
1203 Function   : determinant
1204 Purpose    : Calculate the minor of the given matrix, defined by the columns
1205              specified by values c1, c2, c3
1206 Parameters :
1207 Returns    : determinant value
1208 History    :
1209 ======================================================================== */
1210
1211 static double determinant (const double a[3][4],
1212                            const int    c1,
1213                            const int    c2,
1214                            const int    c3)
1215 {
1216   return a[0][c1]*a[1][c2]*a[2][c3] +
1217          a[0][c2]*a[1][c3]*a[2][c1] +
1218          a[0][c3]*a[1][c1]*a[2][c2] -
1219          a[0][c3]*a[1][c2]*a[2][c1] -
1220          a[0][c2]*a[1][c1]*a[2][c3] -
1221          a[0][c1]*a[1][c3]*a[2][c2];
1222 }
1223
1224 /* =====================================================================
1225 Function   : tri_line_intersect
1226 Purpose    : Intersect a triangle with a line
1227 Parameters : start  - coordinates of the origin of the line
1228              dir    - coordinates of the direction of the line (normalized)
1229              V0     - first vertex of the triangle
1230              V1     - second vertex of the triangle
1231              V2     - third vertex of the triangle
1232              tInter - output value, contains the parameter of the intersection
1233                       point on the line (if found). May be NULL pointer
1234 Returns    : int = 1 if intersection found, 0 otherwise
1235 ======================================================================== */
1236
1237 int NIS_Triangulated::tri_line_intersect (const double      start[3],
1238                                           const double      dir[3],
1239                                           const float       V0[3],
1240                                           const float       V1[3],
1241                                           const float       V2[3],
1242                                           double            * tInter)
1243 {
1244   int res = 0;
1245   const double conf = 1E-15;
1246
1247   const double array[][4] = {
1248     { -dir[0],
1249       double(V1[0] - V0[0]), double(V2[0] - V0[0]), start[0] - double(V0[0]) },
1250     { -dir[1],
1251       double(V1[1] - V0[1]), double(V2[1] - V0[1]), start[1] - double(V0[1]) },
1252     { -dir[2],
1253       double(V1[2] - V0[2]), double(V2[2] - V0[2]), start[2] - double(V0[2]) }
1254   };
1255
1256   const double d  = determinant( array, 0, 1, 2 );
1257   const double dt = determinant( array, 3, 1, 2 );
1258
1259   if (d > conf) {
1260     const double da = determinant( array, 0, 3, 2 );
1261     if (da > -conf) {
1262       const double db = determinant( array, 0, 1, 3 );
1263       res = (db > -conf && da+db <= d+conf);
1264     }
1265   } else if (d < -conf) {
1266     const double da = determinant( array, 0, 3, 2 );
1267     if (da < conf) {
1268       const double db = determinant( array, 0, 1, 3 );
1269       res = (db < conf && da+db >= d-conf);
1270     }
1271   }
1272   if (res && tInter)
1273     *tInter = dt / d;
1274
1275   return res;
1276 }
1277
1278 /* =====================================================================
1279 Function   : tri2d_line_intersect
1280 Purpose    : Intersect a 2D triangle with a 3D line. Z coordinate of triangle
1281            : is zero
1282 Parameters : start  - coordinates of the origin of the line
1283              dir    - coordinates of the direction of the line (normalized)
1284              V0     - first vertex of the triangle
1285              V1     - second vertex of the triangle
1286              V2     - third vertex of the triangle
1287              tInter - output value, contains the parameter of the intersection
1288                       point on the line (if found). May be NULL pointer
1289 Returns    : int = 1 if intersection found, 0 otherwise
1290 ======================================================================== */
1291
1292 int NIS_Triangulated::tri2d_line_intersect (const double      start[3],
1293                                             const double      dir[3],
1294                                             const float       V0[2],
1295                                             const float       V1[2],
1296                                             const float       V2[2],
1297                                             double            * tInter)
1298 {
1299   int res = 0;
1300   const double conf = 1E-15;
1301
1302   // Parallel case is excluded
1303   if (dir[2] * dir[2] > conf)
1304   {
1305     // Find the 2d intersection of (start, dir) with the plane Z = 0.
1306     const double p2d[2] = {
1307       start[0] - dir[0] * start[2] / dir[2],
1308       start[1] - dir[1] * start[2] / dir[2]
1309     };
1310
1311     // Classify the 2d intersection using barycentric coordinates
1312     // (http://www.blackpawn.com/texts/pointinpoly/)
1313     const double v[][2] = {
1314       { static_cast<double>(V1[0]-V0[0]), static_cast<double>(V1[1]-V0[1]) },
1315       { static_cast<double>(V2[0]-V0[0]), static_cast<double>(V2[1]-V0[1]) },
1316       { static_cast<double>(p2d[0]-V0[0]), static_cast<double>(p2d[1]-V0[1]) }
1317     };
1318     const double dot00 = v[0][0]*v[0][0] + v[0][1]*v[0][1];
1319     const double dot01 = v[0][0]*v[1][0] + v[0][1]*v[1][1];
1320     const double dot02 = v[0][0]*v[2][0] + v[0][1]*v[2][1];
1321     const double dot11 = v[1][0]*v[1][0] + v[1][1]*v[1][1];
1322     const double dot12 = v[1][0]*v[2][0] + v[1][1]*v[2][1];
1323     const double denom = (dot00 * dot11 - dot01 * dot01);
1324     if (denom * denom < conf) {
1325       // Point on the 1st side of the triangle
1326       res = (dot01 > -conf && dot00 < dot11 + conf);
1327     } else {
1328       // Barycentric coordinates of the point
1329       const double u = (dot11 * dot02 - dot01 * dot12) / denom;
1330       const double v = (dot00 * dot12 - dot01 * dot02) / denom;
1331       res = (u > -conf) && (v > -conf) && (u + v < 1. + conf);
1332     }
1333   }
1334   if (res && tInter)
1335     *tInter = -start[2] / dir[2];
1336
1337   return res;
1338 }
1339
1340 /* =====================================================================
1341 Function   : seg_line_intersect
1342 Purpose    : Intersect a segment with a line
1343 Parameters : start  - coordinates of the origin of the line
1344              dir    - coordinates of the direction of the line (normalized)
1345              over2  - maximal square distance between the line and the segment
1346                       for intersection state
1347              V0     - first vertex of the segment
1348              V1     - second vertex of the segment
1349              tInter - output value, contains the parameter of the intersection
1350                       point on the line (if found). May be NULL pointer
1351 Returns    : int = 1 if intersection found, 0 otherwise
1352 ======================================================================== */
1353
1354 int NIS_Triangulated::seg_line_intersect (const gp_XYZ&     aStart,
1355                                           const gp_XYZ&     aDir,
1356                                           const double      over2,
1357                                           const float       V0[3],
1358                                           const float       V1[3],
1359                                           double            * tInter)
1360 {
1361   int res = 0;
1362   const gp_XYZ aDirSeg (V1[0]-V0[0], V1[1]-V0[1], V1[2]-V0[2]);
1363   gp_XYZ aDirN = aDirSeg ^ aDir;
1364   Standard_Real aMod2 = aDirN.SquareModulus();
1365   if (aMod2 < Precision::Confusion() * 0.001) {
1366     const gp_XYZ aDelta0 (V0[0]-aStart.X(), V0[1]-aStart.Y(), V0[2]-aStart.Z());
1367     if ((aDelta0 ^ aDir).SquareModulus() < over2) {
1368       res = 1;
1369       const gp_XYZ aDelta1 (V1[0]-aStart.X(),V1[1]-aStart.Y(),V1[2]-aStart.Z());
1370       if (tInter)
1371         * tInter = Min (aDelta0 * aDir, aDelta1 * aDir);
1372     }
1373   } else {
1374     // distance between two unlimited lines
1375     const gp_XYZ aPnt0 (V0[0], V0[1], V0[2]);
1376     const Standard_Real aDistL = (aDirN * aPnt0) - aDirN * aStart;
1377     if (aDistL*aDistL < over2 * aMod2) {
1378       const gp_XYZ aPnt1 (V1[0], V1[1], V1[2]);
1379       Standard_Real aDist[3] = {
1380         ((aPnt0 - aStart) ^ aDir).Modulus(),
1381         ((aPnt1 - aStart) ^ aDir).Modulus(),
1382         0.
1383       };
1384       // Find the intermediate point by interpolation using the distances from
1385       // the end points.
1386       const gp_XYZ aPntI =
1387         (aPnt0 * aDist[1] + aPnt1 * aDist[0]) / (aDist[0] + aDist[1]);
1388       aDist[2] = ((aPntI - aStart) ^ aDir).Modulus();
1389       if (aDist[2] < aDist[0] && aDist[2] < aDist[1]) {
1390         if (aDist[2]*aDist[2] < over2) {
1391           res = 1;
1392           if (tInter)
1393             * tInter = (aPntI - aStart) * aDir;
1394         }
1395       } else if (aDist[0] < aDist[1]) {
1396         if (aDist[0] * aDist[0] < over2) {
1397           res = 1;
1398           if (tInter)
1399             * tInter = (aPnt0 - aStart) * aDir;
1400         }
1401       } else {
1402         if (aDist[1] * aDist[1] < over2) {
1403           res = 1;
1404           if (tInter)
1405             * tInter = (aPnt1 - aStart) * aDir;
1406         }
1407       }
1408     }
1409   }
1410   return res;
1411 }
1412
1413 /* =====================================================================
1414 Function   : seg2d_line_intersect
1415 Purpose    : Intersect a 2D segment with a 3D line
1416 Parameters : start  - coordinates of the origin of the line
1417              dir    - coordinates of the direction of the line (normalized)
1418              over2  - maximal square distance between the line and the segment
1419                       for intersection state
1420              V0     - first vertex of the segment
1421              V1     - second vertex of the segment
1422              tInter - output value, contains the parameter of the intersection
1423                       point on the line (if found). May be NULL pointer
1424 Returns    : int = 1 if intersection found, 0 otherwise
1425 ======================================================================== */
1426
1427 int NIS_Triangulated::seg2d_line_intersect (const gp_XYZ&     aStart,
1428                                             const gp_XYZ&     aDir,
1429                                             const double      over2,
1430                                             const float       V0[2],
1431                                             const float       V1[2],
1432                                             double            * tInter)
1433 {
1434   int res = 0;
1435   const gp_XYZ aDirSeg (V1[0]-V0[0], V1[1]-V0[1], 0.);
1436   gp_XYZ aDirN = aDirSeg ^ aDir;
1437   Standard_Real aMod2 = aDirN.SquareModulus();
1438   if (aMod2 < Precision::Confusion() * 0.001) {
1439     const gp_XYZ aDelta0 (V0[0]-aStart.X(), V0[1]-aStart.Y(), -aStart.Z());
1440     if ((aDelta0 ^ aDir).SquareModulus() < over2) {
1441       res = 1;
1442       const gp_XYZ aDelta1 (V1[0]-aStart.X(), V1[1]-aStart.Y(), -aStart.Z());
1443       if (tInter)
1444         * tInter = Min (aDelta0 * aDir, aDelta1 * aDir);
1445     }
1446   } else {
1447     // distance between two unlimited lines
1448     const gp_XYZ aPnt0 (V0[0], V0[1], 0.);
1449     const Standard_Real aDistL = (aDirN * aPnt0) - aDirN * aStart;
1450     if (aDistL*aDistL < over2 * aMod2) {
1451       const gp_XYZ aPnt1 (V1[0], V1[1], 0.);
1452       Standard_Real aDist[3] = {
1453         ((aPnt0 - aStart) ^ aDir).Modulus(),
1454         ((aPnt1 - aStart) ^ aDir).Modulus(),
1455         0.
1456       };
1457       // Find the intermediate point by interpolation using the distances from
1458       // the end points.
1459       const gp_XYZ aPntI =
1460         (aPnt0 * aDist[1] + aPnt1 * aDist[0]) / (aDist[0] + aDist[1]);
1461       aDist[2] = ((aPntI - aStart) ^ aDir).Modulus();
1462       if (aDist[2] < aDist[0] && aDist[2] < aDist[1]) {
1463         if (aDist[2]*aDist[2] < over2) {
1464           res = 1;
1465           if (tInter)
1466             * tInter = (aPntI - aStart) * aDir;
1467         }
1468       } else if (aDist[0] < aDist[1]) {
1469         if (aDist[0] * aDist[0] < over2) {
1470           res = 1;
1471           if (tInter)
1472             * tInter = (aPnt0 - aStart) * aDir;
1473         }
1474       } else {
1475         if (aDist[1] * aDist[1] < over2) {
1476           res = 1;
1477           if (tInter)
1478             * tInter = (aPnt1 - aStart) * aDir;
1479         }
1480       }
1481     }
1482   }
1483   return res;
1484 }
1485
1486 //=======================================================================
1487 //function : seg_box_intersect
1488 //purpose  :
1489 //=======================================================================
1490
1491 int NIS_Triangulated::seg_box_intersect (const Bnd_B3f& theBox,
1492                                          const gp_Pnt thePnt[2])
1493 {
1494   int aResult (1);
1495   if ((thePnt[1].XYZ() - thePnt[0].XYZ()).SquareModulus()
1496       < Precision::Confusion() * 0.0001)
1497     aResult = 0;
1498   else {
1499     const gp_Dir aDir (thePnt[1].XYZ() - thePnt[0].XYZ());
1500     if (theBox.IsOut (gp_Ax1(thePnt[0], aDir), Standard_True) ||
1501         theBox.IsOut (gp_Ax1(thePnt[1],-aDir), Standard_True))
1502     aResult = 0;
1503   }
1504   return aResult;
1505 }
1506
1507 //=======================================================================
1508 //function : seg_box_included
1509 //purpose  :
1510 //=======================================================================
1511
1512 int NIS_Triangulated::seg_box_included  (const Bnd_B3f& theBox,
1513                                          const gp_Pnt thePnt[2])
1514 {
1515   int aResult (0);
1516   if ((thePnt[1].XYZ() - thePnt[0].XYZ()).SquareModulus()
1517       > Precision::Confusion() * 0.0001)
1518   {
1519     aResult = (theBox.IsOut(thePnt[0].XYZ()) == Standard_False &&
1520                theBox.IsOut(thePnt[1].XYZ()) == Standard_False);
1521   }
1522   return aResult;
1523 }
1524
1525 //=======================================================================
1526 //function : seg_polygon_intersect
1527 //purpose  :
1528 //=======================================================================
1529
1530 int NIS_Triangulated::seg_polygon_intersect
1531                       (const NCollection_List<gp_XY> &thePolygon,
1532                        const gp_XY                    thePnt[2])
1533 {
1534   Standard_Integer aResult = 0;
1535
1536   if (thePolygon.IsEmpty())
1537     return aResult;
1538
1539   gp_XY            aDir(thePnt[1] - thePnt[0]);
1540   Standard_Real    aDist   = aDir.SquareModulus();
1541
1542   if (aDist > aTolConf) {
1543     aDist = Sqrt(aDist);
1544     aDir.Divide(aDist); // Normalize direction.
1545
1546     // Intersect the line passed through thePnt[0] and thePnt[1] with thePolygon
1547     // Line is presented in form Ax + By + C = 0
1548     Standard_Real aA =  aDir.Y();
1549     Standard_Real aB = -aDir.X();
1550     Standard_Real aC = -(aA*thePnt[0].X() + aB*thePnt[0].Y());
1551     gp_XY         aSegment[2];
1552     Standard_Real aSignedD[2];
1553     Standard_Real aDelta = 0.01;
1554
1555     aSegment[0] = thePolygon.Last();
1556     aSignedD[0] = aA*aSegment[0].X() + aB*aSegment[0].Y() + aC;
1557
1558     NCollection_List<gp_XY>::Iterator anIter(thePolygon);
1559
1560     for (; anIter.More(); anIter.Next()) {
1561       aSegment[1] = anIter.Value();
1562       aSignedD[1] = aA*aSegment[1].X() + aB*aSegment[1].Y() + aC;
1563
1564       // Check if there is an intersection.
1565       if (Abs(aSignedD[0]) <= aTolConf || Abs(aSignedD[1]) <= aTolConf) {
1566         Standard_Integer anIndexVtx =
1567                               (Abs(aSignedD[0]) > aTolConf) ? 1 :
1568                               (Abs(aSignedD[1]) > aTolConf) ? 0 : -1;
1569         if (anIndexVtx != -1) {
1570           // Check if the point aSegment[1] is inside the segment.
1571           gp_XY         aDP(aSegment[anIndexVtx] - thePnt[0]);
1572           Standard_Real aParam = aDP.Dot(aDir);
1573
1574           if (aParam >= -aTolConf && aParam <= aDist + aTolConf) {
1575             // Classify a point on the line that is close to aSegment[1] with
1576             // respect to the polygon.
1577             gp_XY aPShift;
1578
1579             if (aParam - aDelta >= 0.) {
1580               aPShift = thePnt[0] + aDir.Multiplied(aParam - aDelta);
1581               if (!IsIn(thePolygon, aPShift))
1582                 aResult = 1;
1583             }
1584
1585             // Try to shift on another direction.
1586             if (!aResult) {
1587               if (aParam + aDelta <= aDist) {
1588                 aPShift = thePnt[0] + aDir.Multiplied(aParam + aDelta);
1589                 if (!IsIn(thePolygon, aPShift))
1590                   aResult = 1;
1591               }
1592             }
1593           }
1594         }
1595       } else if (aSignedD[0]*aSignedD[1] < 0.) {
1596         // Compute intersection of the segment with the line.
1597         gp_XY         aDSeg(aSegment[1] - aSegment[0]);
1598         Standard_Real aSegLen     = aDSeg.Modulus();
1599         Standard_Real aParamOnSeg = aSegLen/(1 + Abs(aSignedD[1]/aSignedD[0]));
1600         gp_XY         aPOnLine
1601                          (aSegment[0] + aDSeg.Multiplied(aParamOnSeg/aSegLen));
1602
1603         // Check if aPOnLine inside the segment thePnt[1] - thePnt[0]
1604         gp_XY         aDP(aPOnLine - thePnt[0]);
1605         Standard_Real aParam = aDP.Dot(aDir);
1606
1607         if (aParam >= -aTolConf && aParam <= aDist + aTolConf) {
1608           gp_XY aPShift;
1609
1610           if (aParam - aDelta >= 0.) {
1611             aPShift = thePnt[0] + aDir.Multiplied(aParam - aDelta);
1612
1613             if (!IsIn(thePolygon, aPShift))
1614               aResult = 1;
1615           }
1616
1617           // Try to shift on another direction.
1618           if (!aResult) {
1619             if (aParam + aDelta <= aDist) {
1620               aPShift = thePnt[0] + aDir.Multiplied(aParam + aDelta);
1621
1622               if (!IsIn(thePolygon, aPShift))
1623                 aResult = 1;
1624             }
1625           }
1626         }
1627       }
1628
1629       if (aResult != 0)
1630         break;
1631
1632       aSegment[0] = aSegment[1];
1633       aSignedD[0] = aSignedD[1];
1634     }
1635   }
1636
1637   return aResult;
1638 }
1639
1640 //=======================================================================
1641 //function : seg_polygon_included
1642 //purpose  :
1643 //=======================================================================
1644
1645 int NIS_Triangulated::seg_polygon_included
1646                       (const NCollection_List<gp_XY> &thePolygon,
1647                        const gp_XY                    thePnt[2])
1648 {
1649   int aResult     = 0;
1650   int anIntersect = seg_polygon_intersect(thePolygon, thePnt);
1651
1652   if (anIntersect == 0) {
1653     if (IsIn(thePolygon, thePnt[0]) && IsIn(thePolygon, thePnt[1]))
1654       aResult = 1;
1655   }
1656
1657   return aResult;
1658 }
1659
1660 //=======================================================================
1661 //function : IsIn
1662 //purpose  :
1663 //=======================================================================
1664
1665 Standard_Boolean NIS_Triangulated::IsIn
1666                       (const NCollection_List<gp_XY> &thePolygon,
1667                        const gp_XY                   &thePoint)
1668 {
1669   if (thePolygon.IsEmpty())
1670     return Standard_False;
1671
1672   Standard_Integer aCounter = 0; // intersections counter
1673   gp_XY            aSegment[2];
1674
1675   aSegment[0] = thePolygon.Last();
1676
1677   NCollection_List<gp_XY>::Iterator anIter(thePolygon);
1678
1679   for (; anIter.More(); anIter.Next()) {
1680     aSegment[1] = anIter.Value();
1681
1682     // Compute projection of the point onto the segment.
1683     Standard_Real       aParam = 0.;
1684     const gp_XY         aDelta = aSegment[1] - aSegment[0];
1685     const gp_XY         aDPP0  = thePoint - aSegment[0];
1686     const Standard_Real aLen2  = aDelta.SquareModulus();
1687
1688     if (IsPositive(aLen2)) {
1689       aParam = (aDelta*aDPP0)/aLen2;
1690
1691       if (aParam < 0.)
1692         aParam = 0.;
1693       else if (aParam > 1.)
1694         aParam = 1.;
1695     }
1696     // Check if the point lies on the segment
1697     gp_XY         aPOnSeg  = aSegment[0]*(1. - aParam) + aSegment[1]*aParam;
1698     Standard_Real aSqrDist = (thePoint - aPOnSeg).SquareModulus();
1699
1700     if (aSqrDist < aTolConf) {
1701       // The point is on the contour.
1702       return Standard_True;
1703     }
1704
1705     // Compute intersection.
1706     const Standard_Real aProd(aDPP0 ^ aDelta);
1707
1708     if (IsPositive(thePoint.X() - aSegment[0].X())) {
1709       if (!IsPositive(thePoint.X() - aSegment[1].X())) {
1710         if (aProd > 0.)
1711           aCounter++;
1712       }
1713     } else {
1714       if (IsPositive(thePoint.X() - aSegment[1].X())) {
1715         if (aProd < 0.)
1716           aCounter++;
1717       }
1718     }
1719
1720     aSegment[0] = aSegment[1];
1721   }
1722
1723   return (aCounter & 0x1);
1724 }
1725
1726 //=======================================================================
1727 //function : allocateNodes
1728 //purpose  :
1729 //=======================================================================
1730
1731 void NIS_Triangulated::allocateNodes (const Standard_Integer nNodes)
1732 {
1733   if (nNodes > 0) {
1734     if (myNNodes > 0)
1735       myAlloc->Free(mypNodes);
1736     myNNodes = nNodes;
1737     mypNodes = static_cast<Standard_ShortReal*>
1738       (myAlloc->Allocate(sizeof(Standard_ShortReal) * myNodeCoord * nNodes));
1739     if (nNodes < 256)
1740       myIndexType = 0;
1741     else if (nNodes < 65536)
1742       myIndexType = 1;
1743     else
1744       myIndexType = 2;
1745   }
1746 }
1747
1748 //=======================================================================
1749 //function : NodeAtInd
1750 //purpose  : Get the node pointed by the i-th index in the array.
1751 //=======================================================================
1752
1753 gp_Pnt NIS_Triangulated::nodeAtInd  (const Standard_Integer * theArray,
1754                                      const Standard_Integer theInd) const
1755 {
1756   if (myIndexType == 0) {
1757     const unsigned char * pInd =
1758       reinterpret_cast<const unsigned char *>(theArray);
1759     return gp_Pnt (mypNodes[myNodeCoord * pInd[theInd] + 0],
1760                    mypNodes[myNodeCoord * pInd[theInd] + 1],
1761                    myNodeCoord < 3 ? 0. :
1762                    mypNodes[myNodeCoord * pInd[theInd] + 2]);
1763   }
1764   if (myIndexType == 1) {
1765     const unsigned short * pInd =
1766       reinterpret_cast<const unsigned short *>(theArray);
1767     return gp_Pnt (mypNodes[myNodeCoord * pInd[theInd] + 0],
1768                    mypNodes[myNodeCoord * pInd[theInd] + 1],
1769                    myNodeCoord < 3 ? 0. :
1770                    mypNodes[myNodeCoord * pInd[theInd] + 2]);
1771   }
1772   return gp_Pnt (mypNodes[myNodeCoord * theArray[theInd] + 0],
1773                  mypNodes[myNodeCoord * theArray[theInd] + 1],
1774                  myNodeCoord < 3 ? 0. :
1775                  mypNodes[myNodeCoord * theArray[theInd] + 2]);
1776 }
1777
1778 //=======================================================================
1779 //function : nodeArrAtInd
1780 //purpose  : Get the node pointed by the i-th index in the array.
1781 //=======================================================================
1782
1783 float* NIS_Triangulated::nodeArrAtInd (const Standard_Integer * theArray,
1784                                        const Standard_Integer theInd) const
1785 {
1786   float* pResult = 0L;
1787   if (myIndexType == 0) {
1788     const unsigned char * pInd =
1789       reinterpret_cast<const unsigned char *>(theArray);
1790     pResult = &mypNodes[myNodeCoord * pInd[theInd]];
1791   }
1792   else if (myIndexType == 1) {
1793     const unsigned short * pInd =
1794       reinterpret_cast<const unsigned short *>(theArray);
1795     pResult = &mypNodes[myNodeCoord * pInd[theInd]];
1796   }
1797   else {
1798     pResult = &mypNodes[myNodeCoord * theArray[theInd]];
1799   }
1800   return pResult;
1801 }