src/VTKViewer/VTKViewer_RectPicker.cxx

   1 //  SALOME VTKViewer : build VTK viewer into Salome desktop
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : VTKViewer_RectPicker.cxx
  25 //  Author : Natalia KOPNOVA
  26 //  Module : SALOME
  27
  28 #include <VTKViewer_RectPicker.h>
  29
  30 #include "vtkActor.h"
  31 #include "vtkAssemblyNode.h"
  32 #include "vtkAssemblyPath.h"
  33 #include "vtkCamera.h"
  34 #include "vtkCommand.h"
  35 #include "vtkImageData.h"
  36 #include "vtkLODProp3D.h"
  37 #include "vtkMapper.h"
  38 #include "vtkMath.h"
  39 #include "vtkObjectFactory.h"
  40 #include "vtkPoints.h"
  41 #include "vtkProp3DCollection.h"
  42 #include "vtkProperty.h"
  43 #include "vtkRenderWindow.h"
  44 #include "vtkRenderer.h"
  45 #include "vtkTransform.h"
  46 #include "vtkVertex.h"
  47 #include "vtkVolume.h"
  48 #include "vtkVolumeMapper.h"
  49
  50 using namespace std;
  51
  52 //----------------------------------------------------------------------------
  53 vtkStandardNewMacro(VTKViewer_RectPicker);
  54 //----------------------------------------------------------------------------
  55
  56 VTKViewer_RectPicker::VTKViewer_RectPicker()
  57 {
  58 }
  59
  60 //----------------------------------------------------------------------------
  61 int VTKViewer_RectPicker::Pick(float selectionX1, float selectionY1, float selectionZ1,
  62                                float selectionX2, float selectionY2, float selectionZ2,
  63                                vtkRenderer *renderer)
  64 {
  65   int k, i;
  66   vtkProp *prop;
  67   vtkCamera *camera;
  68   vtkAbstractMapper3D *mapper = NULL;
  69   float p1World[4][4], p2World[4][4], p1Mapper[4][4], p2Mapper[4][4];
  70   float c1[3], c2[3];
  71   int picked=0;
  72   int *winSize;
  73   float x, y, t, p;
  74   float *viewport;
  75   float cameraPos[4], cameraFP[4];
  76   float *displayCoords, *worldCoords;
  77   float pickPosition[4][3];
  78   double *clipRange;
  79   float ray[4][3], rayLength[4];
  80   int pickable;
  81   int LODId;
  82   float windowLowerLeft[4], windowUpperRight[4];
  83   float bounds[6], tol;
  84   float tF, tB;
  85   float hitPosition[3];
  86   float cameraDOP[3];
  87
  88   //  Initialize picking process
  89   this->Initialize();
  90   this->Renderer = renderer;
  91
  92   /* Selection point is not defined for the rectangle
  93   this->SelectionPoint[0] =
  94   this->SelectionPoint[1] =
  95   this->SelectionPoint[2] =
  96   */
  97
  98   // Invoke start pick method if defined
  99   this->InvokeEvent(vtkCommand::StartPickEvent,NULL);
 100
 101   if ( renderer == NULL )
 102     {
 103     vtkErrorMacro(<<"Must specify renderer!");
 104     return 0;
 105     }
 106
 107   // Get camera focal point and position. Convert to display (screen)
 108   // coordinates. We need a depth value for z-buffer.
 109   //
 110   camera = renderer->GetActiveCamera();
 111   camera->GetPosition((float *)cameraPos); cameraPos[3] = 1.0;
 112   camera->GetFocalPoint((float *)cameraFP); cameraFP[3] = 1.0;
 113
 114   renderer->SetWorldPoint(cameraFP);
 115   renderer->WorldToDisplay();
 116   displayCoords = renderer->GetDisplayPoint();
 117   selectionZ1 = selectionZ2 = displayCoords[2];
 118
 119   // Convert the selection rectangle into world coordinates.
 120   //
 121   renderer->SetDisplayPoint(selectionX1, selectionY1, selectionZ1);
 122   renderer->DisplayToWorld();
 123   worldCoords = renderer->GetWorldPoint();
 124   if ( worldCoords[3] == 0.0 )
 125     {
 126     vtkErrorMacro(<<"Bad homogeneous coordinates");
 127     return 0;
 128     }
 129   for (i=0; i < 3; i++)
 130     {
 131     pickPosition[0][i] = worldCoords[i] / worldCoords[3];
 132     }
 133
 134   renderer->SetDisplayPoint(selectionX1, selectionY2, (selectionZ1+selectionZ2)/2);
 135   renderer->DisplayToWorld();
 136   worldCoords = renderer->GetWorldPoint();
 137   if ( worldCoords[3] == 0.0 )
 138     {
 139     vtkErrorMacro(<<"Bad homogeneous coordinates");
 140     return 0;
 141     }
 142   for (i=0; i < 3; i++)
 143     {
 144     pickPosition[1][i] = worldCoords[i] / worldCoords[3];
 145     }
 146
 147   renderer->SetDisplayPoint(selectionX2, selectionY2, selectionZ2);
 148   renderer->DisplayToWorld();
 149   worldCoords = renderer->GetWorldPoint();
 150   if ( worldCoords[3] == 0.0 )
 151     {
 152     vtkErrorMacro(<<"Bad homogeneous coordinates");
 153     return 0;
 154     }
 155   for (i=0; i < 3; i++)
 156     {
 157     pickPosition[2][i] = worldCoords[i] / worldCoords[3];
 158     }
 159
 160   renderer->SetDisplayPoint(selectionX2, selectionY1, (selectionZ1+selectionZ2)/2);
 161   renderer->DisplayToWorld();
 162   worldCoords = renderer->GetWorldPoint();
 163   if ( worldCoords[3] == 0.0 )
 164     {
 165     vtkErrorMacro(<<"Bad homogeneous coordinates");
 166     return 0;
 167     }
 168   for (i=0; i < 3; i++)
 169     {
 170     pickPosition[3][i] = worldCoords[i] / worldCoords[3];
 171     }
 172
 173   //  Compute the ray endpoints.  The ray is along the line running from
 174   //  the camera position to the selection point, starting where this line
 175   //  intersects the front clipping plane, and terminating where this
 176   //  line intersects the back clipping plane.
 177   for (k=0; k < 4; k++) {
 178     for (i=0; i<3; i++)
 179       {
 180         ray[k][i] = pickPosition[k][i] - cameraPos[i];
 181       }
 182   }
 183   for (i=0; i<3; i++)
 184     {
 185       cameraDOP[i] = cameraFP[i] - cameraPos[i];
 186     }
 187
 188   vtkMath::Normalize(cameraDOP);
 189
 190   for (k=0; k < 4; k++) {
 191     if (( rayLength[k] = vtkMath::Dot(cameraDOP,ray[k])) == 0.0 )
 192       {
 193         vtkWarningMacro("Cannot process points");
 194         return 0;
 195       }
 196   }
 197
 198   clipRange = camera->GetClippingRange();
 199
 200   if ( camera->GetParallelProjection() )
 201     {
 202       for (k=0; k < 4; k++) {
 203         tF = clipRange[0] - rayLength[k];
 204         tB = clipRange[1] - rayLength[k];
 205         for (i=0; i<3; i++)
 206           {
 207             p1World[k][i] = pickPosition[k][i] + tF*cameraDOP[i];
 208             p2World[k][i] = pickPosition[k][i] + tB*cameraDOP[i];
 209           }
 210         p1World[k][3] = p2World[k][3] = 1.0;
 211       }
 212     }
 213   else
 214     {
 215       for (k=0; k < 4; k++) {
 216         tF = clipRange[0] / rayLength[k];
 217         tB = clipRange[1] / rayLength[k];
 218         for (i=0; i<3; i++)
 219           {
 220             p1World[k][i] = cameraPos[i] + tF*ray[k][i];
 221             p2World[k][i] = cameraPos[i] + tB*ray[k][i];
 222           }
 223         p1World[k][3] = p2World[k][3] = 1.0;
 224       }
 225     }
 226
 227   // Compute the center points of ray rectangle
 228   for (i=0; i<3; i++) {
 229     c1[i] = c2[i] = 0;
 230     for (k=0; k<4; k++) {
 231       c1[i] += p1World[k][i];
 232       c2[i] += p2World[k][i];
 233     }
 234     c1[i] = c1[i]/4;
 235     c2[i] = c2[i]/4;
 236   }
 237
 238   // Compute the tolerance in world coordinates.  Do this by
 239   // determining the world coordinates of the diagonal points of the
 240   // window, computing the width of the window in world coordinates, and
 241   // multiplying by the tolerance.
 242   //
 243   viewport = renderer->GetViewport();
 244   winSize = renderer->GetRenderWindow()->GetSize();
 245   x = winSize[0] * viewport[0];
 246   y = winSize[1] * viewport[1];
 247   renderer->SetDisplayPoint(x, y, selectionZ1);
 248   renderer->DisplayToWorld();
 249   renderer->GetWorldPoint(windowLowerLeft);
 250
 251   x = winSize[0] * viewport[2];
 252   y = winSize[1] * viewport[3];
 253   renderer->SetDisplayPoint(x, y, selectionZ2);
 254   renderer->DisplayToWorld();
 255   renderer->GetWorldPoint(windowUpperRight);
 256
 257   for (tol=0.0,i=0; i<3; i++)
 258     {
 259     tol += (windowUpperRight[i] - windowLowerLeft[i]) *
 260       (windowUpperRight[i] - windowLowerLeft[i]);
 261     }
 262
 263   tol = sqrt (tol) * this->Tolerance;
 264
 265   //  Loop over all props.  Transform ray (defined from position of
 266   //  camera to selection point) into coordinates of mapper (not
 267   //  transformed to actors coordinates!  Reduces overall computation!!!).
 268   //  Note that only vtkProp3D's can be picked by vtkPicker.
 269   //
 270   vtkPropCollection *props;
 271   vtkProp *propCandidate;
 272   if ( this->PickFromList )
 273     {
 274     props = this->GetPickList();
 275     }
 276   else
 277     {
 278     props = renderer->GetProps();
 279     }
 280
 281   vtkActor *actor;
 282   vtkLODProp3D *prop3D;
 283   vtkVolume *volume;
 284   vtkAssemblyPath *path;
 285   vtkProperty *tempProperty;
 286   this->Transform->PostMultiply();
 287   for ( props->InitTraversal(); (prop=props->GetNextProp()); )
 288     {
 289     for ( prop->InitPathTraversal(); (path=prop->GetNextPath()); )
 290       {
 291       pickable = 0;
 292       actor = NULL;
 293       propCandidate = path->GetLastNode()->GetProp();
 294       if ( propCandidate->GetPickable() && propCandidate->GetVisibility() )
 295         {
 296         pickable = 1;
 297         if ( (actor=vtkActor::SafeDownCast(propCandidate)) != NULL )
 298           {
 299           mapper = actor->GetMapper();
 300           if ( actor->GetProperty()->GetOpacity() <= 0.0 )
 301             {
 302             pickable = 0;
 303             }
 304           }
 305         else if ( (prop3D=vtkLODProp3D::SafeDownCast(propCandidate)) != NULL )
 306           {
 307           LODId = prop3D->GetPickLODID();
 308           mapper = prop3D->GetLODMapper(LODId);
 309
 310           // if the mapper is a vtkMapper (as opposed to a vtkVolumeMapper),
 311           // then check the transparency to see if the object is pickable
 312           if ( vtkMapper::SafeDownCast(mapper) != NULL)
 313             {
 314             prop3D->GetLODProperty(LODId, &tempProperty);
 315             if ( tempProperty->GetOpacity() <= 0.0 )
 316               {
 317               pickable = 0;
 318               }
 319             }
 320           }
 321         else if ( (volume=vtkVolume::SafeDownCast(propCandidate)) != NULL )
 322           {
 323           mapper = volume->GetMapper();
 324           }
 325         else
 326           {
 327           pickable = 0; //only vtkProp3D's (actors and volumes) can be picked
 328           }
 329         }
 330       //  If actor can be picked, get its composite matrix, invert it, and
 331       //  use the inverted matrix to transform the ray points into mapper
 332       //  coordinates.
 333       if ( pickable  &&  mapper != NULL )
 334         {
 335         vtkMatrix4x4 *LastMatrix = path->GetLastNode()->GetMatrix();
 336         if (LastMatrix == NULL)
 337           {
 338           vtkErrorMacro (<< "Pick: Null matrix.");
 339           return 0;
 340           }
 341         this->Transform->SetMatrix(LastMatrix);
 342         this->Transform->Push();
 343         this->Transform->Inverse();
 344
 345         for (k=0; k < 4; k++) {
 346           this->Transform->TransformPoint(p1World[k],p1Mapper[k]);
 347           this->Transform->TransformPoint(p2World[k],p2Mapper[k]);
 348
 349           for (i=0; i<3; i++)
 350             {
 351               ray[k][i] = p2Mapper[k][i] - p1Mapper[k][i];
 352             }
 353         }
 354
 355         this->Transform->Pop();
 356
 357         //  Have the ray endpoints in mapper space, now need to compare this
 358         //  with the mapper bounds to see whether intersection is possible.
 359         //
 360         //  Get the bounding box of the modeller.  Note that the tolerance is
 361         //  added to the bounding box to make sure things on the edge of the
 362         //  bounding box are picked correctly.
 363         mapper->GetBounds(bounds);
 364         bounds[0] -= tol; bounds[1] += tol;
 365         bounds[2] -= tol; bounds[3] += tol;
 366         bounds[4] -= tol; bounds[5] += tol;
 367         if ( HitBBox(bounds, p1Mapper, ray) ) {
 368           t = this->IntersectWithHex(p1Mapper, p2Mapper, tol, path,
 369                                      (vtkProp3D *)propCandidate, mapper);
 370           if ( t >= 0.0 && t <= 1.0 /*t < VTK_LARGE_FLOAT*/ ) {
 371             picked = 1;
 372             this->Prop3Ds->AddItem((vtkProp3D *)prop);
 373             this->PickedPositions->InsertNextPoint
 374               ((1.0 - t)*c1[0] + t*c2[0],
 375                (1.0 - t)*c1[1] + t*c2[1],
 376                (1.0 - t)*c1[2] + t*c2[2]);
 377
 378             // backwards compatibility: also add to this->Actors
 379             if (actor) {
 380               this->Actors->AddItem(actor);
 381             }
 382           }
 383         }
 384
 385         }//if visible and pickable not transparent and has mapper
 386       }//for all parts
 387     }//for all actors
 388
 389   // Invoke end pick method if defined
 390   this->InvokeEvent(vtkCommand::EndPickEvent,NULL);
 391
 392   return picked;
 393 }
 394
 395 #define SIDE_LEFT 0
 396 #define SIDE_RIGHT 1
 397 #define SIDE_MIDDLE 2
 398
 399 float GetParameterValue(float start, float end, float point)
 400 {
 401   if (start == end) return -VTK_LARGE_FLOAT;
 402   return (point-start)/(end-start);
 403 }
 404
 405 void GetPointCoord(const float start[3], const float end[3], float t, float point[3])
 406 {
 407   int i;
 408   for (i = 0; i < 3; i++) {
 409     point[i] = start[i] + t*(end[i]-start[i]);
 410   }
 411 }
 412
 413 char GetIntersectionPoint(const float start[3], const float end[3],
 414                           const int& index, const float p, float point[3])
 415 {
 416   float t = GetParameterValue(start[index], end[index], p);
 417   char result = 0;
 418   if (t >= 0.0 && t <= 1.0) {
 419     result = 1;
 420     GetPointCoord(start, end, t, point);
 421   }
 422   return result;
 423 }
 424
 425 //----------------------------------------------------------------------------
 426 char VTKViewer_RectPicker::HitBBox (float bounds[6], float origin[4][4], float dir[4][3])
 427 {
 428   int i, j, k, n;
 429   float endray[4][3];
 430
 431   for (k = 0; k < 4; k++) {
 432     for (i = 0; i < 3; i++) {
 433       endray[k][i] = origin[k][i] + dir[k][i];
 434     }
 435   }
 436
 437   // Compute hex bounding box, center point and center direction
 438   float hbounds[6], center[3], ray[3];
 439   for (i = 0; i < 3; i++) {
 440     hbounds[2*i] = hbounds[2*i+1] = origin[0][i];
 441     center[i] = ray[i] = 0;
 442     for (k = 0; k < 4; k++) {
 443       center[i] += origin[k][i];
 444       ray[i] += endray[k][i];
 445       if (origin[k][i] < hbounds[2*i]) {
 446         hbounds[2*i] = origin[k][i];
 447       }
 448       else if (origin[k][i] > hbounds[2*i+1])
 449         hbounds[2*i+1] = origin[k][i];
 450       if (endray[k][i] < hbounds[2*i])
 451         hbounds[2*i] = endray[k][i];
 452       else if (endray[k][i] > hbounds[2*i+1])
 453         hbounds[2*i+1] = endray[k][i];
 454     }
 455     center[i] = center[i]/4;
 456     ray[i] = ray[i]/4;
 457     ray[i] = ray[i] - center[i];
 458   }
 459
 460   // Check for intersection between bouning boxes
 461   for (i = 0; i < 3; i++) {
 462     if (bounds[2*i+1] < hbounds[2*i] || bounds[2*i] > hbounds[2*i+1])
 463       return 0;
 464   }
 465
 466   // Check if one of the origin point lays inside bbox
 467   char inside;
 468   for (k = 0; k < 4; k++) {
 469     inside = 1;
 470     for (i = 0; i < 3; i++) {
 471       if (origin[k][i] < bounds[2*i] || origin[k][i] > bounds[2*i+1]) {
 472         inside = 0;
 473         break;
 474       }
 475     }
 476     if (inside) return 1;
 477   }
 478
 479   // Find the closest coord plane for the center point
 480   char side[3];
 481   float coordPlane[3];
 482   inside = 1;
 483   for (i = 0; i < 3; i++) {
 484     if (center[i] < bounds[2*i]) {
 485       inside = 0;
 486       coordPlane[i] = bounds[2*i];
 487       side[i] = SIDE_LEFT;
 488     }
 489     else if (center[i] > bounds[2*i+1]) {
 490       inside = 0;
 491       coordPlane[i] = bounds[2*i+1];
 492       side[i] = SIDE_RIGHT;
 493     }
 494     else {
 495       coordPlane[i] = (ray[i]<0.0) ? bounds[2*i] : bounds[2*i+1];
 496       side[i] = SIDE_MIDDLE;
 497     }
 498   }
 499   if (inside) return 1;
 500
 501   // Calculate parametric distances to the planes and find the max
 502   float maxT[3];
 503   int whichPlane = 0;
 504   char defined = 0;
 505   for (i = 0; i < 3; i++) {
 506     if (side[i] != SIDE_MIDDLE && ray[i] != 0.0) {
 507       maxT[i] = (coordPlane[i]-center[i])/ray[i];
 508       defined = 1;
 509     }
 510     else
 511       maxT[i] = -1.0;
 512   }
 513   for (i = 0; i < 3; i++) {
 514     if (maxT[whichPlane] < maxT[i])
 515       whichPlane = i;
 516   }
 517
 518   // Check for intersection along the center ray
 519   float coord;
 520   if (maxT[whichPlane] <= 1.0 && maxT[whichPlane] >= 0.0) {
 521     inside = 1;
 522     for (i = 0; i < 3; i++) {
 523       if (i != whichPlane) {
 524         coord = center[i] + maxT[whichPlane]*ray[i];
 525         if (coord < bounds[2*i] || coord > bounds[2*i+1])
 526           inside = 0;
 527       }
 528     }
 529     if (inside) return 1;
 530   }
 531
 532   // Define the intersection plane
 533   if (!defined) {
 534     for (i = 0; i < 3; i++) {
 535       if (ray[i] != 0.0) {
 536         maxT[i] = (coordPlane[i]-center[i])/ray[i];
 537       }
 538       else
 539         maxT[i] = VTK_LARGE_FLOAT;
 540     }
 541     for (i = 0; i < 3; i++) {
 542       if (maxT[whichPlane] > maxT[i])
 543         whichPlane = i;
 544     }
 545   }
 546
 547   // Compute the intersection between hex and coord plane
 548   float t[4];
 549   for (k = 0; k < 4; k++) {
 550     if (dir[k][whichPlane] != 0.0) {
 551       t[k] = (coordPlane[whichPlane]-origin[k][whichPlane])/dir[k][whichPlane];
 552     }
 553     else {
 554       t[k] = VTK_LARGE_FLOAT;
 555     }
 556   }
 557
 558   vtkPoints* aPoints = vtkPoints::New();
 559   float p[3], q[3], t1;
 560   for (k = 0; k < 4; k++) {
 561     n = (k+1)%4; // next point
 562     if (t[k] > 1.0) {
 563       if (t[n] < 1.0) {
 564         // find intersection point
 565         t1 = GetParameterValue(endray[k][whichPlane], endray[n][whichPlane], coordPlane[whichPlane]);
 566         if (t1 > 0.0 && t1 < 1.0) {
 567           GetPointCoord(endray[k], endray[n], t1, p);
 568           aPoints->InsertNextPoint(p[0], p[1], p[2]);
 569         }
 570       }
 571       if (t[n] < 0.0) {
 572         // find second intersection point
 573         t1 = GetParameterValue(origin[k][whichPlane], origin[n][whichPlane], coordPlane[whichPlane]);
 574         if (t1 > 0.0 && t1 < 1.0) {
 575           GetPointCoord(origin[k], origin[n], t1, p);
 576           aPoints->InsertNextPoint(p[0], p[1], p[2]);
 577         }
 578       }
 579     }
 580     else if (t[k] < 0.0) {
 581       if (t[n] > 0.0) {
 582         // find intersection point
 583         t1 = GetParameterValue(origin[k][whichPlane], origin[n][whichPlane], coordPlane[whichPlane]);
 584         if (t1 > 0.0 && t1 < 1.0) {
 585           GetPointCoord(origin[k], origin[n], t1, p);
 586           aPoints->InsertNextPoint(p[0], p[1], p[2]);
 587         }
 588       }
 589     }
 590     else {
 591       // find intersection point
 592       GetPointCoord(origin[k], endray[k], t[k], p);
 593       aPoints->InsertNextPoint(p[0], p[1], p[2]);
 594
 595       if (t[n] < 0.0) {
 596         // find second intersection point
 597         t1 = GetParameterValue(origin[k][whichPlane], origin[n][whichPlane], coordPlane[whichPlane]);
 598         if (t1 > 0.0 && t1 < 1.0) {
 599           GetPointCoord(origin[k], origin[n], t1, p);
 600           aPoints->InsertNextPoint(p[0], p[1], p[2]);
 601         }
 602       }
 603       else if (t[n] > 1.0) {
 604         // find second intersection point
 605         t1 = GetParameterValue(endray[k][whichPlane], endray[n][whichPlane], coordPlane[whichPlane]);
 606         if (t1 > 0.0 && t1 < 1.0) {
 607           GetPointCoord(endray[k], endray[n], t1, p);
 608           aPoints->InsertNextPoint(p[0], p[1], p[2]);
 609         }
 610       }
 611     }
 612   }
 613   n = aPoints->GetNumberOfPoints();
 614   if (n == 0) {
 615     aPoints->Delete();
 616     return 0;
 617   }
 618
 619   if (n == 1) {
 620     aPoints->GetPoint(0, p);
 621     inside = 1;
 622     for (i = 0; i < 3; i++) {
 623       if (i != whichPlane) {
 624         if (p[i] < bounds[2*i] || p[i] > bounds[2*i+1]) {
 625           inside = 0; break;
 626         }
 627       }
 628     }
 629     aPoints->Delete();
 630     return inside;
 631   }
 632
 633   // Analize intersection
 634   int nearPlane, boundPlane = -1;
 635   float boundCoord, boundMin, boundMax;
 636   char intersect = 0;
 637   for (k = 0; k < n; k++) {
 638     aPoints->GetPoint(k, p);
 639     j = k+1; if (j == n) j = 0;
 640     aPoints->GetPoint(j, q);
 641     inside = 1;
 642     nearPlane = 0;
 643     // if the point is inside bbox
 644     for (i = 0; i < 3; i++) {
 645       if (i != whichPlane) {
 646         if (p[i] < bounds[2*i]) {
 647           side[i] = SIDE_LEFT;
 648           maxT[i] = GetParameterValue(p[i], q[i], bounds[2*i]);
 649           inside = 0;
 650         }
 651         else if (p[i] > bounds[2*i+1]) {
 652           side[i] = SIDE_RIGHT;
 653           maxT[i] = GetParameterValue(p[i], q[i], bounds[2*i+1]);
 654           inside = 0;
 655         }
 656         else {
 657           side[i] = SIDE_MIDDLE;
 658           maxT[i] = -1.0;
 659         }
 660       }
 661       else maxT[i] = -1.0;
 662       if (maxT[i] > maxT[nearPlane]) nearPlane = i;
 663     }
 664     if (inside) break;
 665     // if segment intersects bbox
 666     if (maxT[nearPlane] >= 0.0 && maxT[nearPlane] <= 1.0) {
 667       for (i = 0; i < 3; i++) {
 668         if (i != whichPlane && i != nearPlane) {
 669           coord = p[i] + maxT[nearPlane]*(q[i]-p[i]);
 670           if (coord >= bounds[2*i] && coord <= bounds[2*i+1]) {
 671             intersect = 1; break;
 672           }
 673         }
 674       }
 675       // intersect with boundPlane
 676       if (boundPlane == -1) {
 677         boundCoord = p[nearPlane] + maxT[nearPlane]*(q[nearPlane]-p[nearPlane]);
 678         boundPlane = nearPlane;
 679         for (i = 0; i < 3; i++) {
 680           if (i != whichPlane && i != boundPlane) {
 681             coord = p[i] + maxT[nearPlane]*(q[i]-p[i]);
 682             boundMin = boundMax = coord;
 683           }
 684         }
 685       }
 686       else {
 687         t1 = GetParameterValue(p[boundPlane], q[boundPlane], boundCoord);
 688         if (t1 >= 0.0 && t1 <= 1.0) {
 689           for (i = 0; i < 3; i++) {
 690             if (i != whichPlane && i != boundPlane) {
 691               coord = p[i] + t1*(q[i]-p[i]);
 692               if (coord < boundMin) boundMin = coord;
 693               if (coord > boundMax) boundMax = coord;
 694             }
 695           }
 696         }
 697       }
 698     }
 699     if (intersect) break;
 700   }
 701   aPoints->Delete();
 702   if (inside || intersect) {
 703     return 1;
 704   }
 705
 706   inside = 1;
 707   for (i = 0; i < 3; i++) {
 708     if (i != whichPlane && i != boundPlane) {
 709       if (boundMin > bounds[2*i+1] || boundMax < bounds[2*i])
 710         inside = 0;
 711     }
 712   }
 713
 714   return inside;
 715 }
 716
 717 //----------------------------------------------------------------------------
 718 char VTKViewer_RectPicker::PointInside (float p[3], float p1[4][4], float p2[4][4], float tol)
 719 {
 720   int i, j, k;
 721   float t, coord[3];
 722
 723   // Fix one coordinate (x, for example) and
 724   // compute intersection with coordinate plane
 725   vtkPoints* aPoints = vtkPoints::New();
 726   int mode = 0;
 727   for (k = 0; k < 4; k++) {
 728     j = k+1; if (j == 4) j = 0;
 729     switch (mode) {
 730     case 0:
 731       if (GetIntersectionPoint(p1[k], p1[j], 0, p[0], coord)) {
 732         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 733         mode = 0;
 734       }
 735       if (GetIntersectionPoint(p1[k], p2[k], 0, p[0], coord)) {
 736         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 737         mode = 1;
 738       }
 739       if (GetIntersectionPoint(p2[k], p2[j], 0, p[0], coord)) {
 740         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 741         mode = 2;
 742       }
 743       /*
 744       if ((p1[k][0]-p[0])*(p2[k][0]-p[0]) <= 0) {
 745         t = GetParameterValue(p1[k][0], p2[k][0], p[0]);
 746         if (t >= 0.0 && t <= 1.0) {
 747           GetPointCoord(p1[k], p2[k], t, coord);
 748           aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 749         }
 750       }
 751       */
 752       break;
 753     case 1:
 754       if (GetIntersectionPoint(p1[k], p2[k], 0, p[0], coord)) {
 755         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 756         mode = 1;
 757       }
 758       if (GetIntersectionPoint(p2[k], p2[j], 0, p[0], coord)) {
 759         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 760         mode = 2;
 761       }
 762       if (GetIntersectionPoint(p1[k], p1[j], 0, p[0], coord)) {
 763         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 764         mode = 0;
 765       }
 766       /*
 767       if ((p1[k][0]-p[0])*(p1[j][0]-p[0]) <= 0) {
 768         t = GetParameterValue(p1[k][0], p1[j][0], p[0]);
 769         if (t > 0.0 && t < 1.0) {
 770           GetPointCoord(p1[k], p1[j], t, coord);
 771           aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 772         }
 773       }
 774       */
 775       break;
 776     case 2:
 777       if (GetIntersectionPoint(p2[k], p2[j], 0, p[0], coord)) {
 778         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 779         mode = 2;
 780       }
 781       if (GetIntersectionPoint(p1[k], p2[k], 0, p[0], coord)) {
 782         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 783         mode = 1;
 784       }
 785       if (GetIntersectionPoint(p1[k], p1[j], 0, p[0], coord)) {
 786         aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 787         mode = 0;
 788       }
 789       /*
 790       if ((p2[k][0]-p[0])*(p2[j][0]-p[0]) <= 0) {
 791         t = GetParameterValue(p2[k][0], p2[j][0], p[0]);
 792         if (t > 0.0 && t < 1.0) {
 793           GetPointCoord(p2[k], p2[j], t, coord);
 794           aPoints->InsertNextPoint(coord[0], coord[1], coord[2]);
 795         }
 796       }
 797       */
 798       break;
 799     }
 800   }
 801   int n = aPoints->GetNumberOfPoints();
 802   //cout << "---> Points in X projection " << n << endl;
 803   if (n == 0) {
 804     aPoints->Delete();
 805     return 0;
 806   }
 807
 808   // Fix the second coord and define bounds
 809   float zMin = VTK_LARGE_FLOAT, zMax = -VTK_LARGE_FLOAT, z, ncoord[3];
 810   char inside = 0;
 811   for (k = 0; k < n; k++) {
 812     aPoints->GetPoint(k, coord);
 813     //cout << "  P" << k << " (" << coord[0] << ", " << coord[1] << ", " << coord[2] << ")";
 814     j = k+1; if (j == n) j = 0;
 815     if (j == k) {
 816       if (p[1] == coord[1] && p[2] == coord[2]) {
 817         inside = 1;
 818       }
 819       break;
 820     }
 821     aPoints->GetPoint(j, ncoord);
 822     t = GetParameterValue(coord[1], ncoord[1], p[1]);
 823     if (t >= 0.0 && t <= 1) {
 824       z = coord[2] + t*(ncoord[2]-coord[2]);
 825       if (z < zMin) zMin = z;
 826       if (z > zMax) zMax = z;
 827     }
 828   }
 829   //cout << endl << " Zmin = " << zMin << ", Zmax = " << zMax  << endl;
 830   if (!inside) {
 831     if (p[2] <= (zMax+tol) && p[2] >= (zMin-tol))
 832       inside = 1;
 833   }
 834
 835   aPoints->Delete();
 836   return inside;
 837 }
 838
 839 //----------------------------------------------------------------------------
 840 float VTKViewer_RectPicker::IntersectWithHex(float p1[4][4], float p2[4][4], float tol,
 841                                              vtkAssemblyPath *path, vtkProp3D *prop3D,
 842                                              vtkAbstractMapper3D *mapper)
 843 {
 844   int i, k;
 845   float *center, p0[3], ray[3], rayFactor, t;
 846
 847   // Get the data from the modeler
 848   //
 849   center = mapper->GetCenter();
 850
 851   if (!PointInside(center, p1, p2)) {
 852     return 2.0;
 853   }
 854
 855   //   Determine appropriate info
 856   //
 857   for (i = 0; i < 3; i++) {
 858     p0[i] = ray[i] = 0;
 859     for (k = 0; k < 4; k++) {
 860       p0[i] += p1[k][i];
 861       ray[i] += p2[k][i];
 862     }
 863     p0[i] = p0[i]/4;
 864     ray[i] = ray[i]/4;
 865     ray[i] = ray[i] - p0[i];
 866   }
 867   if (( rayFactor = vtkMath::Dot(ray,ray)) == 0.0 ) {
 868     vtkErrorMacro("Cannot process points");
 869     return 2.0;
 870   }
 871
 872   // Project the center point onto the ray and determine its parametric value
 873   //
 874   t = (ray[0]*(center[0]-p0[0]) + ray[1]*(center[1]-p0[1])
 875        + ray[2]*(center[2]-p0[2])) / rayFactor;
 876
 877   if ( t >= 0.0 && t <= 1.0 && t < this->GlobalTMin ) {
 878     this->MarkPicked(path, prop3D, mapper, t, center);
 879   }
 880   return t;
 881 }