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 under
   8 // the terms of the GNU Lesser General Public License 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 nbSteps = (unsigned)myNNodes * myNodeCoord;
 776     for (unsigned int iNode = 0; iNode < nbSteps; iNode += myNodeCoord)
 777     {
 778       gp_XYZ aPnt (mypNodes[iNode+0], mypNodes[iNode+1], 0.);
 779       if (myNodeCoord > 2)
 780         aPnt.SetZ (mypNodes[iNode+2]);
 781       theTrf.Transforms (aPnt);
 782       if (theBox.IsOut (aPnt)) {
 783         if (isFullIn) {
 784           aResult = Standard_False;
 785           break;
 786         }
 787       } else {
 788         if (isFullIn == Standard_False) {
 789           aResult = Standard_True;
 790           break;
 791         }
 792       }
 793     }
 794   }
 795   if (aResult == isFullIn) {
 796     if (myType & Type_Segments) {
 797       for (Standard_Integer i = 0; i < myNLineNodes; i+=2) {
 798         const gp_Pnt aPnt[2] = {
 799           nodeAtInd(mypLines, i+0).Transformed(theTrf),
 800           nodeAtInd(mypLines, i+1).Transformed(theTrf)
 801         };
 802         if (isFullIn) {
 803           if (seg_box_included (theBox, aPnt) == 0) {
 804             aResult = Standard_False;
 805             break;
 806           }
 807         } else {
 808           if (seg_box_intersect (theBox, aPnt)) {
 809             aResult = Standard_True;
 810             break;
 811           }
 812         }
 813       }
 814     } else if (myType & Type_Line) {
 815       for (Standard_Integer i = 1; i < myNLineNodes; i++) {
 816         const gp_Pnt aPnt[2] = {
 817           nodeAtInd(mypLines, i-1).Transformed(theTrf),
 818           nodeAtInd(mypLines, i+0).Transformed(theTrf)
 819         };
 820         if (isFullIn) {
 821           if (seg_box_included (theBox, aPnt) == 0) {
 822             aResult = Standard_False;
 823             break;
 824           }
 825         } else {
 826           if (seg_box_intersect (theBox, aPnt)) {
 827             aResult = Standard_True;
 828             break;
 829           }
 830         }
 831       }
 832       if (aResult == isFullIn && (myType & Type_Loop)) {
 833         const gp_Pnt aPntLast[2] = {
 834           nodeAtInd(mypLines, myNLineNodes-1).Transformed(theTrf),
 835           nodeAtInd(mypLines, 0).Transformed(theTrf)
 836         };
 837         if (isFullIn) {
 838           if (seg_box_included (theBox, aPntLast) == 0)
 839             aResult = Standard_False;
 840         } else {
 841           if (seg_box_intersect (theBox, aPntLast))
 842             aResult = Standard_True;
 843         }
 844       }
 845     } else if ((myType & Type_Polygons) && myIsDrawPolygons) {
 846       for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
 847         const Standard_Integer * aPoly = mypPolygons[iPoly];
 848         const Standard_Integer nNodes = NPolygonNodes(iPoly);
 849         for (Standard_Integer i = 1; i < nNodes; i++) {
 850           // index is incremented by 1 for the head number in the array
 851           const gp_Pnt aPnt[2] = {
 852             nodeAtInd(aPoly, i+0).Transformed(theTrf),
 853             nodeAtInd(aPoly, i+1).Transformed(theTrf)
 854           };
 855           if (isFullIn) {
 856             if (seg_box_included (theBox, aPnt) == 0) {
 857               aResult = Standard_False;
 858               break;
 859             }
 860           } else {
 861             if (seg_box_intersect (theBox, aPnt)) {
 862               aResult = Standard_True;
 863               break;
 864             }
 865           }
 866         }
 867         if (aResult == isFullIn) {
 868           const gp_Pnt aPntLast[2] = {
 869             nodeAtInd(aPoly, nNodes).Transformed(theTrf),
 870             nodeAtInd(aPoly, 1).Transformed(theTrf)
 871           };
 872           if (isFullIn) {
 873             if (seg_box_included (theBox, aPntLast) == 0)
 874               aResult = Standard_False;
 875           } else {
 876             if (seg_box_intersect (theBox, aPntLast))
 877               aResult = Standard_True;
 878           }
 879         }
 880       }
 881     }
 882   }
 883   return aResult;
 884 }
 885
 886 //=======================================================================
 887 //function : Intersect
 888 //purpose  :
 889 //=======================================================================
 890
 891 Standard_Real NIS_Triangulated::Intersect (const gp_Ax1&       theAxis,
 892                                            const Standard_Real theOver) const
 893 {
 894   Standard_Real aResult (RealLast());
 895   Standard_Real start[3], dir[3];
 896   theAxis.Location().Coord(start[0], start[1], start[2]);
 897   theAxis.Direction().Coord(dir[0], dir[1], dir[2]);
 898   double anInter;
 899
 900   if ((myType & Type_Triangulation) && (myIsDrawPolygons == Standard_False))
 901     for (Standard_Integer i = 0; i < myNTriangles; i++) {
 902       Standard_Boolean isIntersect(Standard_False);
 903       if (myIndexType == 0) {
 904         const unsigned char * pTri =
 905           reinterpret_cast<unsigned char *>(mypTriangles) + (3 * i);
 906         if (myNodeCoord > 2)
 907           isIntersect = tri_line_intersect (start, dir,
 908                                             &mypNodes[myNodeCoord * pTri[0]],
 909                                             &mypNodes[myNodeCoord * pTri[1]],
 910                                             &mypNodes[myNodeCoord * pTri[2]],
 911                                             &anInter);
 912         else
 913           isIntersect = tri2d_line_intersect (start, dir,
 914                                               &mypNodes[myNodeCoord * pTri[0]],
 915                                               &mypNodes[myNodeCoord * pTri[1]],
 916                                               &mypNodes[myNodeCoord * pTri[2]],
 917                                               &anInter);
 918       } else if (myIndexType == 1) {
 919         const unsigned short * pTri =
 920           reinterpret_cast<unsigned short *>(mypTriangles) + (3 * i);
 921         if (myNodeCoord > 2)
 922           isIntersect = tri_line_intersect (start, dir,
 923                                             &mypNodes[myNodeCoord * pTri[0]],
 924                                             &mypNodes[myNodeCoord * pTri[1]],
 925                                             &mypNodes[myNodeCoord * pTri[2]],
 926                                             &anInter);
 927         else
 928           isIntersect = tri2d_line_intersect (start, dir,
 929                                               &mypNodes[myNodeCoord * pTri[0]],
 930                                               &mypNodes[myNodeCoord * pTri[1]],
 931                                               &mypNodes[myNodeCoord * pTri[2]],
 932                                               &anInter);
 933       } else {
 934         const Standard_Integer * pTri = &mypTriangles[3 * i];
 935         if (myNodeCoord > 2)
 936           isIntersect = tri_line_intersect (start, dir,
 937                                             &mypNodes[myNodeCoord * pTri[0]],
 938                                             &mypNodes[myNodeCoord * pTri[1]],
 939                                             &mypNodes[myNodeCoord * pTri[2]],
 940                                             &anInter);
 941         else
 942           isIntersect = tri2d_line_intersect (start, dir,
 943                                               &mypNodes[myNodeCoord * pTri[0]],
 944                                               &mypNodes[myNodeCoord * pTri[1]],
 945                                               &mypNodes[myNodeCoord * pTri[2]],
 946                                               &anInter);
 947       }
 948       if (isIntersect && anInter < aResult)
 949         aResult = anInter;
 950     }
 951
 952   const Standard_Real anOver2 = theOver*theOver;
 953   if (myType & Type_Segments) {
 954     Standard_Integer i = 0;
 955     if (myNodeCoord > 2)
 956       for (; i < myNLineNodes; i+=2) {
 957         if (seg_line_intersect (theAxis.Location().XYZ(),
 958                                 theAxis.Direction().XYZ(), anOver2,
 959                                 nodeArrAtInd(mypLines, i+0),
 960                                 nodeArrAtInd(mypLines, i+1),
 961                                 &anInter))
 962           if (anInter < aResult)
 963             aResult = anInter;
 964       }
 965     else
 966       for (; i < myNLineNodes; i+=2) {
 967         if (seg2d_line_intersect (theAxis.Location().XYZ(),
 968                                   theAxis.Direction().XYZ(), anOver2,
 969                                   nodeArrAtInd(mypLines, i+0),
 970                                   nodeArrAtInd(mypLines, i+1),
 971                                   &anInter))
 972           if (anInter < aResult)
 973             aResult = anInter;
 974       }
 975   } else if (myType & Type_Line) {
 976     Standard_Integer i = 1;
 977     if (myNodeCoord > 2) {
 978       for (; i < myNLineNodes; i++) {
 979         if (seg_line_intersect (theAxis.Location().XYZ(),
 980                                 theAxis.Direction().XYZ(), anOver2,
 981                                 nodeArrAtInd(mypLines, i-1),
 982                                 nodeArrAtInd(mypLines, i-0),
 983                                 &anInter))
 984           if (anInter < aResult)
 985             aResult = anInter;
 986       }
 987       if (myType & Type_Loop)
 988         if (seg_line_intersect (theAxis.Location().XYZ(),
 989                                 theAxis.Direction().XYZ(), anOver2,
 990                                 nodeArrAtInd(mypLines, myNLineNodes-1),
 991                                 nodeArrAtInd(mypLines, 0),
 992                                 &anInter))
 993           if (anInter < aResult)
 994             aResult = anInter;
 995     } else {
 996       for (; i < myNLineNodes; i++) {
 997         if (seg2d_line_intersect (theAxis.Location().XYZ(),
 998                                   theAxis.Direction().XYZ(), anOver2,
 999                                   nodeArrAtInd(mypLines, i-1),
1000                                   nodeArrAtInd(mypLines, i-0),
1001                                   &anInter))
1002           if (anInter < aResult)
1003             aResult = anInter;
1004       }
1005       if (myType & Type_Loop)
1006         if (seg2d_line_intersect (theAxis.Location().XYZ(),
1007                                   theAxis.Direction().XYZ(), anOver2,
1008                                   nodeArrAtInd(mypLines, myNLineNodes-1),
1009                                   nodeArrAtInd(mypLines, 0),
1010                                   &anInter))
1011           if (anInter < aResult)
1012             aResult = anInter;
1013     }
1014   }
1015
1016   if ((myType & Type_Polygons) && myIsDrawPolygons) {
1017     for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
1018       const Standard_Integer * aPoly = mypPolygons[iPoly];
1019       const Standard_Integer nNodes = NPolygonNodes(iPoly);
1020       Standard_Integer i = 1;
1021       if (myNodeCoord > 2) {
1022         for (; i < nNodes; i++) {
1023           // Node index is incremented for the head of polygon indice array
1024           if (seg_line_intersect (theAxis.Location().XYZ(),
1025                                   theAxis.Direction().XYZ(), anOver2,
1026                                   nodeArrAtInd(aPoly, i+0),
1027                                   nodeArrAtInd(aPoly, i+1),
1028                                   &anInter))
1029             if (anInter < aResult)
1030               aResult = anInter;
1031         }
1032         if (seg_line_intersect (theAxis.Location().XYZ(),
1033                                 theAxis.Direction().XYZ(), anOver2,
1034                                 nodeArrAtInd(aPoly, nNodes),
1035                                 nodeArrAtInd(aPoly, 1),
1036                                 &anInter))
1037           if (anInter < aResult)
1038             aResult = anInter;
1039       } else {
1040         for (; i < nNodes; i++) {
1041           // Node index is incremented for the head of polygon indice array
1042           if (seg2d_line_intersect (theAxis.Location().XYZ(),
1043                                     theAxis.Direction().XYZ(), anOver2,
1044                                     nodeArrAtInd(aPoly, i+0),
1045                                     nodeArrAtInd(aPoly, i+1),
1046                                     &anInter))
1047             if (anInter < aResult)
1048               aResult = anInter;
1049         }
1050         if (seg2d_line_intersect (theAxis.Location().XYZ(),
1051                                   theAxis.Direction().XYZ(), anOver2,
1052                                   nodeArrAtInd(aPoly, nNodes),
1053                                   nodeArrAtInd(aPoly, 1),
1054                                   &anInter))
1055           if (anInter < aResult)
1056             aResult = anInter;
1057       }
1058     }
1059   }
1060   return aResult;
1061 }
1062
1063 //=======================================================================
1064 //function : Intersect
1065 //purpose  :
1066 //=======================================================================
1067 Standard_Boolean NIS_Triangulated::Intersect
1068                     (const NCollection_List<gp_XY> &thePolygon,
1069                      const gp_Trsf                 &theTrf,
1070                      const Standard_Boolean         isFullIn) const
1071 {
1072   Standard_Boolean aResult (isFullIn);
1073
1074   if ((myType & Type_Triangulation) && myIsDrawPolygons == Standard_False) {
1075     unsigned int nbSteps = (unsigned)myNNodes * myNodeCoord;
1076     for (unsigned int iNode = 0; iNode < nbSteps; iNode += myNodeCoord)
1077     {
1078       gp_XYZ aPnt (mypNodes[iNode+0], mypNodes[iNode+1], 0.);
1079       if (myNodeCoord > 2)
1080         aPnt.SetZ (mypNodes[iNode+2]);
1081       theTrf.Transforms (aPnt);
1082
1083       gp_XY aP2d(aPnt.X(), aPnt.Y());
1084
1085       if (!NIS_Triangulated::IsIn(thePolygon, aP2d)) {
1086         if (isFullIn) {
1087           aResult = Standard_False;
1088           break;
1089         }
1090       } else {
1091         if (isFullIn == Standard_False) {
1092           aResult = Standard_True;
1093           break;
1094         }
1095       }
1096     }
1097   }
1098   if (aResult == isFullIn) {
1099     if (myType & Type_Segments) {
1100       for (Standard_Integer i = 0; i < myNLineNodes; i+=2) {
1101         const gp_Pnt aPnt[2] = {
1102           nodeAtInd(mypLines, i+0).Transformed(theTrf),
1103           nodeAtInd(mypLines, i+1).Transformed(theTrf)
1104         };
1105         const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1106                                 gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1107
1108         if (isFullIn) {
1109           if (seg_polygon_included (thePolygon, aP2d) == 0) {
1110             aResult = Standard_False;
1111             break;
1112           }
1113         } else {
1114           if (seg_polygon_intersect (thePolygon, aP2d)) {
1115             aResult = Standard_True;
1116             break;
1117           }
1118         }
1119       }
1120     } else if (myType & Type_Line) {
1121       for (Standard_Integer i = 1; i < myNLineNodes; i++) {
1122         const gp_Pnt aPnt[2] = {
1123           nodeAtInd(mypLines, i-1).Transformed(theTrf),
1124           nodeAtInd(mypLines, i+0).Transformed(theTrf)
1125         };
1126         const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1127                                 gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1128         if (isFullIn) {
1129           if (seg_polygon_included (thePolygon, aP2d) == 0) {
1130             aResult = Standard_False;
1131             break;
1132           }
1133         } else {
1134           if (seg_polygon_intersect (thePolygon, aP2d)) {
1135             aResult = Standard_True;
1136             break;
1137           }
1138         }
1139       }
1140       if (aResult == isFullIn && (myType & Type_Loop)) {
1141         const gp_Pnt aPntLast[2] = {
1142           nodeAtInd(mypLines, myNLineNodes-1).Transformed(theTrf),
1143           nodeAtInd(mypLines, 0).Transformed(theTrf)
1144         };
1145         const gp_XY aP2dLast[2] = { gp_XY(aPntLast[0].X(), aPntLast[0].Y()),
1146                                     gp_XY(aPntLast[1].X(), aPntLast[1].Y()) };
1147
1148         if (isFullIn) {
1149           if (seg_polygon_included (thePolygon, aP2dLast) == 0)
1150             aResult = Standard_False;
1151         } else {
1152           if (seg_polygon_intersect (thePolygon, aP2dLast))
1153             aResult = Standard_True;
1154         }
1155       }
1156     } else if ((myType & Type_Polygons) && myIsDrawPolygons) {
1157       for (unsigned int iPoly = 0; iPoly < myNPolygons; iPoly++) {
1158         const Standard_Integer * aPoly = mypPolygons[iPoly];
1159         const Standard_Integer nNodes = NPolygonNodes(iPoly);
1160         for (Standard_Integer i = 1; i < nNodes; i++) {
1161           const gp_Pnt aPnt[2] = {
1162             nodeAtInd(aPoly, i+0).Transformed(theTrf),
1163             nodeAtInd(aPoly, i+1).Transformed(theTrf)
1164           };
1165           const gp_XY aP2d[2] = { gp_XY(aPnt[0].X(), aPnt[0].Y()),
1166                                   gp_XY(aPnt[1].X(), aPnt[1].Y()) };
1167           if (isFullIn) {
1168             if (seg_polygon_included (thePolygon, aP2d) == 0) {
1169               aResult = Standard_False;
1170               break;
1171             }
1172           } else {
1173             if (seg_polygon_intersect (thePolygon, aP2d)) {
1174               aResult = Standard_True;
1175               break;
1176             }
1177           }
1178         }
1179         if (aResult == isFullIn) {
1180           const gp_Pnt aPntLast[2] = {
1181             nodeAtInd(aPoly, nNodes).Transformed(theTrf),
1182             nodeAtInd(aPoly, 1).Transformed(theTrf)
1183           };
1184           const gp_XY aP2dLast[2] = { gp_XY(aPntLast[0].X(), aPntLast[0].Y()),
1185                                       gp_XY(aPntLast[1].X(), aPntLast[1].Y()) };
1186           if (isFullIn) {
1187             if (seg_polygon_included (thePolygon, aP2dLast) == 0)
1188               aResult = Standard_False;
1189           } else {
1190             if (seg_polygon_intersect (thePolygon, aP2dLast))
1191               aResult = Standard_True;
1192           }
1193         }
1194       }
1195     }
1196   }
1197   return aResult;
1198 }
1199
1200 /* =====================================================================
1201 Function   : determinant
1202 Purpose    : Calculate the minor of the given matrix, defined by the columns
1203              specified by values c1, c2, c3
1204 Parameters :
1205 Returns    : determinant value
1206 History    :
1207 ======================================================================== */
1208
1209 static double determinant (const double a[3][4],
1210                            const int    c1,
1211                            const int    c2,
1212                            const int    c3)
1213 {
1214   return a[0][c1]*a[1][c2]*a[2][c3] +
1215          a[0][c2]*a[1][c3]*a[2][c1] +
1216          a[0][c3]*a[1][c1]*a[2][c2] -
1217          a[0][c3]*a[1][c2]*a[2][c1] -
1218          a[0][c2]*a[1][c1]*a[2][c3] -
1219          a[0][c1]*a[1][c3]*a[2][c2];
1220 }
1221
1222 /* =====================================================================
1223 Function   : tri_line_intersect
1224 Purpose    : Intersect a triangle with a line
1225 Parameters : start  - coordinates of the origin of the line
1226              dir    - coordinates of the direction of the line (normalized)
1227              V0     - first vertex of the triangle
1228              V1     - second vertex of the triangle
1229              V2     - third vertex of the triangle
1230              tInter - output value, contains the parameter of the intersection
1231                       point on the line (if found). May be NULL pointer
1232 Returns    : int = 1 if intersection found, 0 otherwise
1233 ======================================================================== */
1234
1235 int NIS_Triangulated::tri_line_intersect (const double      start[3],
1236                                           const double      dir[3],
1237                                           const float       V0[3],
1238                                           const float       V1[3],
1239                                           const float       V2[3],
1240                                           double            * tInter)
1241 {
1242   int res = 0;
1243   const double conf = 1E-15;
1244
1245   const double array[][4] = {
1246     { -dir[0],
1247       double(V1[0] - V0[0]), double(V2[0] - V0[0]), start[0] - double(V0[0]) },
1248     { -dir[1],
1249       double(V1[1] - V0[1]), double(V2[1] - V0[1]), start[1] - double(V0[1]) },
1250     { -dir[2],
1251       double(V1[2] - V0[2]), double(V2[2] - V0[2]), start[2] - double(V0[2]) }
1252   };
1253
1254   const double d  = determinant( array, 0, 1, 2 );
1255   const double dt = determinant( array, 3, 1, 2 );
1256
1257   if (d > conf) {
1258     const double da = determinant( array, 0, 3, 2 );
1259     if (da > -conf) {
1260       const double db = determinant( array, 0, 1, 3 );
1261       res = (db > -conf && da+db <= d+conf);
1262     }
1263   } else if (d < -conf) {
1264     const double da = determinant( array, 0, 3, 2 );
1265     if (da < conf) {
1266       const double db = determinant( array, 0, 1, 3 );
1267       res = (db < conf && da+db >= d-conf);
1268     }
1269   }
1270   if (res && tInter)
1271     *tInter = dt / d;
1272
1273   return res;
1274 }
1275
1276 /* =====================================================================
1277 Function   : tri2d_line_intersect
1278 Purpose    : Intersect a 2D triangle with a 3D line. Z coordinate of triangle
1279            : is zero
1280 Parameters : start  - coordinates of the origin of the line
1281              dir    - coordinates of the direction of the line (normalized)
1282              V0     - first vertex of the triangle
1283              V1     - second vertex of the triangle
1284              V2     - third vertex of the triangle
1285              tInter - output value, contains the parameter of the intersection
1286                       point on the line (if found). May be NULL pointer
1287 Returns    : int = 1 if intersection found, 0 otherwise
1288 ======================================================================== */
1289
1290 int NIS_Triangulated::tri2d_line_intersect (const double      start[3],
1291                                             const double      dir[3],
1292                                             const float       V0[2],
1293                                             const float       V1[2],
1294                                             const float       V2[2],
1295                                             double            * tInter)
1296 {
1297   int res = 0;
1298   const double conf = 1E-15;
1299
1300   // Parallel case is excluded
1301   if (dir[2] * dir[2] > conf)
1302   {
1303     // Find the 2d intersection of (start, dir) with the plane Z = 0.
1304     const double p2d[2] = {
1305       start[0] - dir[0] * start[2] / dir[2],
1306       start[1] - dir[1] * start[2] / dir[2]
1307     };
1308
1309     // Classify the 2d intersection using barycentric coordinates
1310     // (http://www.blackpawn.com/texts/pointinpoly/)
1311     const double v[][2] = {
1312       { static_cast<double>(V1[0]-V0[0]), static_cast<double>(V1[1]-V0[1]) },
1313       { static_cast<double>(V2[0]-V0[0]), static_cast<double>(V2[1]-V0[1]) },
1314       { static_cast<double>(p2d[0]-V0[0]), static_cast<double>(p2d[1]-V0[1]) }
1315     };
1316     const double dot00 = v[0][0]*v[0][0] + v[0][1]*v[0][1];
1317     const double dot01 = v[0][0]*v[1][0] + v[0][1]*v[1][1];
1318     const double dot02 = v[0][0]*v[2][0] + v[0][1]*v[2][1];
1319     const double dot11 = v[1][0]*v[1][0] + v[1][1]*v[1][1];
1320     const double dot12 = v[1][0]*v[2][0] + v[1][1]*v[2][1];
1321     const double denom = (dot00 * dot11 - dot01 * dot01);
1322     if (denom * denom < conf) {
1323       // Point on the 1st side of the triangle
1324       res = (dot01 > -conf && dot00 < dot11 + conf);
1325     } else {
1326       // Barycentric coordinates of the point
1327       const double u = (dot11 * dot02 - dot01 * dot12) / denom;
1328       const double v = (dot00 * dot12 - dot01 * dot02) / denom;
1329       res = (u > -conf) && (v > -conf) && (u + v < 1. + conf);
1330     }
1331   }
1332   if (res && tInter)
1333     *tInter = -start[2] / dir[2];
1334
1335   return res;
1336 }
1337
1338 /* =====================================================================
1339 Function   : seg_line_intersect
1340 Purpose    : Intersect a segment with a line
1341 Parameters : start  - coordinates of the origin of the line
1342              dir    - coordinates of the direction of the line (normalized)
1343              over2  - maximal square distance between the line and the segment
1344                       for intersection state
1345              V0     - first vertex of the segment
1346              V1     - second vertex of the segment
1347              tInter - output value, contains the parameter of the intersection
1348                       point on the line (if found). May be NULL pointer
1349 Returns    : int = 1 if intersection found, 0 otherwise
1350 ======================================================================== */
1351
1352 int NIS_Triangulated::seg_line_intersect (const gp_XYZ&     aStart,
1353                                           const gp_XYZ&     aDir,
1354                                           const double      over2,
1355                                           const float       V0[3],
1356                                           const float       V1[3],
1357                                           double            * tInter)
1358 {
1359   int res = 0;
1360   const gp_XYZ aDirSeg (V1[0]-V0[0], V1[1]-V0[1], V1[2]-V0[2]);
1361   gp_XYZ aDirN = aDirSeg ^ aDir;
1362   Standard_Real aMod2 = aDirN.SquareModulus();
1363   if (aMod2 < Precision::Confusion() * 0.001) {
1364     const gp_XYZ aDelta0 (V0[0]-aStart.X(), V0[1]-aStart.Y(), V0[2]-aStart.Z());
1365     if ((aDelta0 ^ aDir).SquareModulus() < over2) {
1366       res = 1;
1367       const gp_XYZ aDelta1 (V1[0]-aStart.X(),V1[1]-aStart.Y(),V1[2]-aStart.Z());
1368       if (tInter)
1369         * tInter = Min (aDelta0 * aDir, aDelta1 * aDir);
1370     }
1371   } else {
1372     // distance between two unlimited lines
1373     const gp_XYZ aPnt0 (V0[0], V0[1], V0[2]);
1374     const Standard_Real aDistL = (aDirN * aPnt0) - aDirN * aStart;
1375     if (aDistL*aDistL < over2 * aMod2) {
1376       const gp_XYZ aPnt1 (V1[0], V1[1], V1[2]);
1377       Standard_Real aDist[3] = {
1378         ((aPnt0 - aStart) ^ aDir).Modulus(),
1379         ((aPnt1 - aStart) ^ aDir).Modulus(),
1380         0.
1381       };
1382       // Find the intermediate point by interpolation using the distances from
1383       // the end points.
1384       const gp_XYZ aPntI =
1385         (aPnt0 * aDist[1] + aPnt1 * aDist[0]) / (aDist[0] + aDist[1]);
1386       aDist[2] = ((aPntI - aStart) ^ aDir).Modulus();
1387       if (aDist[2] < aDist[0] && aDist[2] < aDist[1]) {
1388         if (aDist[2]*aDist[2] < over2) {
1389           res = 1;
1390           if (tInter)
1391             * tInter = (aPntI - aStart) * aDir;
1392         }
1393       } else if (aDist[0] < aDist[1]) {
1394         if (aDist[0] * aDist[0] < over2) {
1395           res = 1;
1396           if (tInter)
1397             * tInter = (aPnt0 - aStart) * aDir;
1398         }
1399       } else {
1400         if (aDist[1] * aDist[1] < over2) {
1401           res = 1;
1402           if (tInter)
1403             * tInter = (aPnt1 - aStart) * aDir;
1404         }
1405       }
1406     }
1407   }
1408   return res;
1409 }
1410
1411 /* =====================================================================
1412 Function   : seg2d_line_intersect
1413 Purpose    : Intersect a 2D segment with a 3D line
1414 Parameters : start  - coordinates of the origin of the line
1415              dir    - coordinates of the direction of the line (normalized)
1416              over2  - maximal square distance between the line and the segment
1417                       for intersection state
1418              V0     - first vertex of the segment
1419              V1     - second vertex of the segment
1420              tInter - output value, contains the parameter of the intersection
1421                       point on the line (if found). May be NULL pointer
1422 Returns    : int = 1 if intersection found, 0 otherwise
1423 ======================================================================== */
1424
1425 int NIS_Triangulated::seg2d_line_intersect (const gp_XYZ&     aStart,
1426                                             const gp_XYZ&     aDir,
1427                                             const double      over2,
1428                                             const float       V0[2],
1429                                             const float       V1[2],
1430                                             double            * tInter)
1431 {
1432   int res = 0;
1433   const gp_XYZ aDirSeg (V1[0]-V0[0], V1[1]-V0[1], 0.);
1434   gp_XYZ aDirN = aDirSeg ^ aDir;
1435   Standard_Real aMod2 = aDirN.SquareModulus();
1436   if (aMod2 < Precision::Confusion() * 0.001) {
1437     const gp_XYZ aDelta0 (V0[0]-aStart.X(), V0[1]-aStart.Y(), -aStart.Z());
1438     if ((aDelta0 ^ aDir).SquareModulus() < over2) {
1439       res = 1;
1440       const gp_XYZ aDelta1 (V1[0]-aStart.X(), V1[1]-aStart.Y(), -aStart.Z());
1441       if (tInter)
1442         * tInter = Min (aDelta0 * aDir, aDelta1 * aDir);
1443     }
1444   } else {
1445     // distance between two unlimited lines
1446     const gp_XYZ aPnt0 (V0[0], V0[1], 0.);
1447     const Standard_Real aDistL = (aDirN * aPnt0) - aDirN * aStart;
1448     if (aDistL*aDistL < over2 * aMod2) {
1449       const gp_XYZ aPnt1 (V1[0], V1[1], 0.);
1450       Standard_Real aDist[3] = {
1451         ((aPnt0 - aStart) ^ aDir).Modulus(),
1452         ((aPnt1 - aStart) ^ aDir).Modulus(),
1453         0.
1454       };
1455       // Find the intermediate point by interpolation using the distances from
1456       // the end points.
1457       const gp_XYZ aPntI =
1458         (aPnt0 * aDist[1] + aPnt1 * aDist[0]) / (aDist[0] + aDist[1]);
1459       aDist[2] = ((aPntI - aStart) ^ aDir).Modulus();
1460       if (aDist[2] < aDist[0] && aDist[2] < aDist[1]) {
1461         if (aDist[2]*aDist[2] < over2) {
1462           res = 1;
1463           if (tInter)
1464             * tInter = (aPntI - aStart) * aDir;
1465         }
1466       } else if (aDist[0] < aDist[1]) {
1467         if (aDist[0] * aDist[0] < over2) {
1468           res = 1;
1469           if (tInter)
1470             * tInter = (aPnt0 - aStart) * aDir;
1471         }
1472       } else {
1473         if (aDist[1] * aDist[1] < over2) {
1474           res = 1;
1475           if (tInter)
1476             * tInter = (aPnt1 - aStart) * aDir;
1477         }
1478       }
1479     }
1480   }
1481   return res;
1482 }
1483
1484 //=======================================================================
1485 //function : seg_box_intersect
1486 //purpose  :
1487 //=======================================================================
1488
1489 int NIS_Triangulated::seg_box_intersect (const Bnd_B3f& theBox,
1490                                          const gp_Pnt thePnt[2])
1491 {
1492   int aResult (1);
1493   if ((thePnt[1].XYZ() - thePnt[0].XYZ()).SquareModulus()
1494       < Precision::Confusion() * 0.0001)
1495     aResult = 0;
1496   else {
1497     const gp_Dir aDir (thePnt[1].XYZ() - thePnt[0].XYZ());
1498     if (theBox.IsOut (gp_Ax1(thePnt[0], aDir), Standard_True) ||
1499         theBox.IsOut (gp_Ax1(thePnt[1],-aDir), Standard_True))
1500     aResult = 0;
1501   }
1502   return aResult;
1503 }
1504
1505 //=======================================================================
1506 //function : seg_box_included
1507 //purpose  :
1508 //=======================================================================
1509
1510 int NIS_Triangulated::seg_box_included  (const Bnd_B3f& theBox,
1511                                          const gp_Pnt thePnt[2])
1512 {
1513   int aResult (0);
1514   if ((thePnt[1].XYZ() - thePnt[0].XYZ()).SquareModulus()
1515       > Precision::Confusion() * 0.0001)
1516   {
1517     aResult = (theBox.IsOut(thePnt[0].XYZ()) == Standard_False &&
1518                theBox.IsOut(thePnt[1].XYZ()) == Standard_False);
1519   }
1520   return aResult;
1521 }
1522
1523 //=======================================================================
1524 //function : seg_polygon_intersect
1525 //purpose  :
1526 //=======================================================================
1527
1528 int NIS_Triangulated::seg_polygon_intersect
1529                       (const NCollection_List<gp_XY> &thePolygon,
1530                        const gp_XY                    thePnt[2])
1531 {
1532   Standard_Integer aResult = 0;
1533
1534   if (thePolygon.IsEmpty())
1535     return aResult;
1536
1537   gp_XY            aDir(thePnt[1] - thePnt[0]);
1538   Standard_Real    aDist   = aDir.SquareModulus();
1539
1540   if (aDist > aTolConf) {
1541     aDist = Sqrt(aDist);
1542     aDir.Divide(aDist); // Normalize direction.
1543
1544     // Intersect the line passed through thePnt[0] and thePnt[1] with thePolygon
1545     // Line is presented in form Ax + By + C = 0
1546     Standard_Real aA =  aDir.Y();
1547     Standard_Real aB = -aDir.X();
1548     Standard_Real aC = -(aA*thePnt[0].X() + aB*thePnt[0].Y());
1549     gp_XY         aSegment[2];
1550     Standard_Real aSignedD[2];
1551     Standard_Real aDelta = 0.01;
1552
1553     aSegment[0] = thePolygon.Last();
1554     aSignedD[0] = aA*aSegment[0].X() + aB*aSegment[0].Y() + aC;
1555
1556     NCollection_List<gp_XY>::Iterator anIter(thePolygon);
1557
1558     for (; anIter.More(); anIter.Next()) {
1559       aSegment[1] = anIter.Value();
1560       aSignedD[1] = aA*aSegment[1].X() + aB*aSegment[1].Y() + aC;
1561
1562       // Check if there is an intersection.
1563       if (Abs(aSignedD[0]) <= aTolConf || Abs(aSignedD[1]) <= aTolConf) {
1564         Standard_Integer anIndexVtx =
1565                               (Abs(aSignedD[0]) > aTolConf) ? 1 :
1566                               (Abs(aSignedD[1]) > aTolConf) ? 0 : -1;
1567         if (anIndexVtx != -1) {
1568           // Check if the point aSegment[1] is inside the segment.
1569           gp_XY         aDP(aSegment[anIndexVtx] - thePnt[0]);
1570           Standard_Real aParam = aDP.Dot(aDir);
1571
1572           if (aParam >= -aTolConf && aParam <= aDist + aTolConf) {
1573             // Classify a point on the line that is close to aSegment[1] with
1574             // respect to the polygon.
1575             gp_XY aPShift;
1576
1577             if (aParam - aDelta >= 0.) {
1578               aPShift = thePnt[0] + aDir.Multiplied(aParam - aDelta);
1579               if (!IsIn(thePolygon, aPShift))
1580                 aResult = 1;
1581             }
1582
1583             // Try to shift on another direction.
1584             if (!aResult) {
1585               if (aParam + aDelta <= aDist) {
1586                 aPShift = thePnt[0] + aDir.Multiplied(aParam + aDelta);
1587                 if (!IsIn(thePolygon, aPShift))
1588                   aResult = 1;
1589               }
1590             }
1591           }
1592         }
1593       } else if (aSignedD[0]*aSignedD[1] < 0.) {
1594         // Compute intersection of the segment with the line.
1595         gp_XY         aDSeg(aSegment[1] - aSegment[0]);
1596         Standard_Real aSegLen     = aDSeg.Modulus();
1597         Standard_Real aParamOnSeg = aSegLen/(1 + Abs(aSignedD[1]/aSignedD[0]));
1598         gp_XY         aPOnLine
1599                          (aSegment[0] + aDSeg.Multiplied(aParamOnSeg/aSegLen));
1600
1601         // Check if aPOnLine inside the segment thePnt[1] - thePnt[0]
1602         gp_XY         aDP(aPOnLine - thePnt[0]);
1603         Standard_Real aParam = aDP.Dot(aDir);
1604
1605         if (aParam >= -aTolConf && aParam <= aDist + aTolConf) {
1606           gp_XY aPShift;
1607
1608           if (aParam - aDelta >= 0.) {
1609             aPShift = thePnt[0] + aDir.Multiplied(aParam - aDelta);
1610
1611             if (!IsIn(thePolygon, aPShift))
1612               aResult = 1;
1613           }
1614
1615           // Try to shift on another direction.
1616           if (!aResult) {
1617             if (aParam + aDelta <= aDist) {
1618               aPShift = thePnt[0] + aDir.Multiplied(aParam + aDelta);
1619
1620               if (!IsIn(thePolygon, aPShift))
1621                 aResult = 1;
1622             }
1623           }
1624         }
1625       }
1626
1627       if (aResult != 0)
1628         break;
1629
1630       aSegment[0] = aSegment[1];
1631       aSignedD[0] = aSignedD[1];
1632     }
1633   }
1634
1635   return aResult;
1636 }
1637
1638 //=======================================================================
1639 //function : seg_polygon_included
1640 //purpose  :
1641 //=======================================================================
1642
1643 int NIS_Triangulated::seg_polygon_included
1644                       (const NCollection_List<gp_XY> &thePolygon,
1645                        const gp_XY                    thePnt[2])
1646 {
1647   int aResult     = 0;
1648   int anIntersect = seg_polygon_intersect(thePolygon, thePnt);
1649
1650   if (anIntersect == 0) {
1651     if (IsIn(thePolygon, thePnt[0]) && IsIn(thePolygon, thePnt[1]))
1652       aResult = 1;
1653   }
1654
1655   return aResult;
1656 }
1657
1658 //=======================================================================
1659 //function : IsIn
1660 //purpose  :
1661 //=======================================================================
1662
1663 Standard_Boolean NIS_Triangulated::IsIn
1664                       (const NCollection_List<gp_XY> &thePolygon,
1665                        const gp_XY                   &thePoint)
1666 {
1667   if (thePolygon.IsEmpty())
1668     return Standard_False;
1669
1670   Standard_Integer aCounter = 0; // intersections counter
1671   gp_XY            aSegment[2];
1672
1673   aSegment[0] = thePolygon.Last();
1674
1675   NCollection_List<gp_XY>::Iterator anIter(thePolygon);
1676
1677   for (; anIter.More(); anIter.Next()) {
1678     aSegment[1] = anIter.Value();
1679
1680     // Compute projection of the point onto the segment.
1681     Standard_Real       aParam = 0.;
1682     const gp_XY         aDelta = aSegment[1] - aSegment[0];
1683     const gp_XY         aDPP0  = thePoint - aSegment[0];
1684     const Standard_Real aLen2  = aDelta.SquareModulus();
1685
1686     if (IsPositive(aLen2)) {
1687       aParam = (aDelta*aDPP0)/aLen2;
1688
1689       if (aParam < 0.)
1690         aParam = 0.;
1691       else if (aParam > 1.)
1692         aParam = 1.;
1693     }
1694     // Check if the point lies on the segment
1695     gp_XY         aPOnSeg  = aSegment[0]*(1. - aParam) + aSegment[1]*aParam;
1696     Standard_Real aSqrDist = (thePoint - aPOnSeg).SquareModulus();
1697
1698     if (aSqrDist < aTolConf) {
1699       // The point is on the contour.
1700       return Standard_True;
1701     }
1702
1703     // Compute intersection.
1704     const Standard_Real aProd(aDPP0 ^ aDelta);
1705
1706     if (IsPositive(thePoint.X() - aSegment[0].X())) {
1707       if (!IsPositive(thePoint.X() - aSegment[1].X())) {
1708         if (aProd > 0.)
1709           aCounter++;
1710       }
1711     } else {
1712       if (IsPositive(thePoint.X() - aSegment[1].X())) {
1713         if (aProd < 0.)
1714           aCounter++;
1715       }
1716     }
1717
1718     aSegment[0] = aSegment[1];
1719   }
1720
1721   return (aCounter & 0x1);
1722 }
1723
1724 //=======================================================================
1725 //function : allocateNodes
1726 //purpose  :
1727 //=======================================================================
1728
1729 void NIS_Triangulated::allocateNodes (const Standard_Integer nNodes)
1730 {
1731   if (nNodes > 0) {
1732     if (myNNodes > 0)
1733       myAlloc->Free(mypNodes);
1734     myNNodes = nNodes;
1735     mypNodes = static_cast<Standard_ShortReal*>
1736       (myAlloc->Allocate(sizeof(Standard_ShortReal) * myNodeCoord * nNodes));
1737     if (nNodes < 256)
1738       myIndexType = 0;
1739     else if (nNodes < 65536)
1740       myIndexType = 1;
1741     else
1742       myIndexType = 2;
1743   }
1744 }
1745
1746 //=======================================================================
1747 //function : NodeAtInd
1748 //purpose  : Get the node pointed by the i-th index in the array.
1749 //=======================================================================
1750
1751 gp_Pnt NIS_Triangulated::nodeAtInd  (const Standard_Integer * theArray,
1752                                      const Standard_Integer theInd) const
1753 {
1754   if (myIndexType == 0) {
1755     const unsigned char * pInd =
1756       reinterpret_cast<const unsigned char *>(theArray);
1757     return gp_Pnt (mypNodes[myNodeCoord * pInd[theInd] + 0],
1758                    mypNodes[myNodeCoord * pInd[theInd] + 1],
1759                    myNodeCoord < 3 ? 0. :
1760                    mypNodes[myNodeCoord * pInd[theInd] + 2]);
1761   }
1762   if (myIndexType == 1) {
1763     const unsigned short * pInd =
1764       reinterpret_cast<const unsigned short *>(theArray);
1765     return gp_Pnt (mypNodes[myNodeCoord * pInd[theInd] + 0],
1766                    mypNodes[myNodeCoord * pInd[theInd] + 1],
1767                    myNodeCoord < 3 ? 0. :
1768                    mypNodes[myNodeCoord * pInd[theInd] + 2]);
1769   }
1770   return gp_Pnt (mypNodes[myNodeCoord * theArray[theInd] + 0],
1771                  mypNodes[myNodeCoord * theArray[theInd] + 1],
1772                  myNodeCoord < 3 ? 0. :
1773                  mypNodes[myNodeCoord * theArray[theInd] + 2]);
1774 }
1775
1776 //=======================================================================
1777 //function : nodeArrAtInd
1778 //purpose  : Get the node pointed by the i-th index in the array.
1779 //=======================================================================
1780
1781 float* NIS_Triangulated::nodeArrAtInd (const Standard_Integer * theArray,
1782                                        const Standard_Integer theInd) const
1783 {
1784   float* pResult = 0L;
1785   if (myIndexType == 0) {
1786     const unsigned char * pInd =
1787       reinterpret_cast<const unsigned char *>(theArray);
1788     pResult = &mypNodes[myNodeCoord * pInd[theInd]];
1789   }
1790   else if (myIndexType == 1) {
1791     const unsigned short * pInd =
1792       reinterpret_cast<const unsigned short *>(theArray);
1793     pResult = &mypNodes[myNodeCoord * pInd[theInd]];
1794   }
1795   else {
1796     pResult = &mypNodes[myNodeCoord * theArray[theInd]];
1797   }
1798   return pResult;
1799 }