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