src/VTKViewer/VTKViewer_CellRectPicker.cxx

   1 #include <VTKViewer_CellRectPicker.h>
   2
   3 #include "vtkGenericCell.h"
   4 #include "vtkAssemblyNode.h"
   5 #include "vtkAssemblyPath.h"
   6 #include "vtkObjectFactory.h"
   7 #include "vtkImageData.h"
   8 #include "vtkMapper.h"
   9 #include "vtkVolumeMapper.h"
  10 #include "vtkMath.h"
  11 #include "vtkLine.h"
  12 #include "vtkQuad.h"
  13 #include "vtkPoints.h"
  14 #include "vtkPlane.h"
  15 #include "vtkActor.h"
  16
  17 vtkStandardNewMacro(VTKViewer_CellRectPicker);
  18
  19 /*!Constructor.\n
  20  * \li create instance of generic cell
  21  * \li create two instances of Quad.
  22  */
  23 VTKViewer_CellRectPicker::VTKViewer_CellRectPicker()
  24 {
  25   this->Cell = vtkGenericCell::New();
  26   this->Quad1 = vtkQuad::New();
  27   this->Quad2 = vtkQuad::New();
  28 }
  29
  30 /*!Destructor.
  31  * \li Delete generic cell
  32  * \li Delete two quads.
  33  */
  34 VTKViewer_CellRectPicker::~VTKViewer_CellRectPicker()
  35 {
  36   this->Cell->Delete();
  37   this->Quad1->Delete();
  38   this->Quad2->Delete();
  39 }
  40
  41 /*!Clear \a ActorData and call VTKViewer_RectPicker::Initialize() method*/
  42 void VTKViewer_CellRectPicker::Initialize()
  43 {
  44   ActorData.clear();
  45   this->VTKViewer_RectPicker::Initialize();
  46 }
  47
  48 /*!*/
  49 float VTKViewer_CellRectPicker::IntersectWithHex(float p1[4][4], float p2[4][4], float tol,
  50                                                  vtkAssemblyPath *path, vtkProp3D *prop3D,
  51                                                  vtkAbstractMapper3D *m)
  52 {
  53   vtkDataSet *input;
  54   vtkMapper *mapper;
  55   vtkVolumeMapper *volumeMapper;
  56
  57   // Get the underlying dataset
  58   if ( (mapper=vtkMapper::SafeDownCast(m)) != NULL ) {
  59     input = mapper->GetInput();
  60   }
  61   else if ( (volumeMapper=vtkVolumeMapper::SafeDownCast(m)) != NULL ) {
  62     input = volumeMapper->GetInput();
  63   }
  64   else {
  65     return VTK_LARGE_FLOAT;
  66   }
  67
  68   vtkIdType numCells;
  69   if ( (numCells = input->GetNumberOfCells()) < 1 ) {
  70     return 2.0;
  71   }
  72
  73   int i, k, subId;
  74   float bounds[6], center[3], cp1[3], cp2[3];
  75   float pcoord[3], weight[4], dist;
  76   float t1[2], t2[2], x[3], t;
  77
  78   for (k = 0; k < 4; k++) {
  79     this->Quad1->PointIds->SetId(k, k);
  80     this->Quad1->Points->SetPoint(k, p1[k]);
  81     this->Quad2->PointIds->SetId(k, k);
  82     this->Quad2->Points->SetPoint(k, p2[k]);
  83   }
  84
  85   /*
  86   cout << "---> Selection area:"   << endl;
  87   for (k = 0; k < 4; k++)
  88   cout << "\t(" << p1[k][0] << ", " << p1[k][1] << ", " << p1[k][2] << ")";
  89   cout << endl;
  90   for (k = 0; k < 4; k++)
  91   cout << "\t(" << p2[k][0] << ", " << p2[k][1] << ", " << p2[k][2] << ")";
  92   cout << endl;
  93   */
  94
  95   vtkIdType cellId;
  96   VTKViewer_CellData cellData;
  97   VTKViewer_CellDataSet dataList;
  98
  99   char inside; int n; float *point;
 100   float tMin = VTK_LARGE_FLOAT, xMin[3];
 101   for (cellId = 0; cellId < numCells; cellId++) {
 102     input->GetCell(cellId, this->Cell);
 103
 104     this->Cell->GetBounds(bounds);
 105     for (i = 0; i < 3; i++)
 106       center[i] = (bounds[2*i]+bounds[2*i+1])/2;
 107
 108       /*
 109       if (!VTKViewer_RectPicker::PointInside(center, p1, p2, 0.01))
 110       continue;
 111     */
 112
 113     inside = 1;
 114     n = this->Cell->Points->GetNumberOfPoints();
 115     //cout << "---> Info cell " << cellId << " (" << n << "): " << endl;
 116     for (k = 0; k < n; k++) {
 117       point = this->Cell->Points->GetPoint(k);
 118       //cout << "  P (" << point[0] << ", " << point[1] << ", " << point[2] << ")";
 119       if (!VTKViewer_RectPicker::PointInside(point, p1, p2)) {
 120         inside = 0; break;
 121       }
 122     }
 123     //cout << endl;
 124     if (!inside) continue;
 125
 126     //cout << "---> Inside cell " << cellId << endl;
 127
 128     this->Quad1->EvaluatePosition(center, 0, subId, pcoord, dist, weight);
 129     t1[0] = pcoord[0]; t1[1] = pcoord[1];
 130     this->Quad2->EvaluatePosition(center, 0, subId, pcoord, dist, weight);
 131     t2[0] = pcoord[0]; t2[1] = pcoord[1];
 132
 133     pcoord[0] = (t1[0]+t2[0])/2; pcoord[1] = (t1[1]+t2[1])/2; pcoord[2] = 0;
 134
 135     this->Quad1->EvaluateLocation(subId, pcoord, cp1, weight);
 136     this->Quad2->EvaluateLocation(subId, pcoord, cp2, weight);
 137
 138     if (this->Cell->IntersectWithLine(cp1, cp2, tol, t, x, pcoord, subId)) {
 139       cellData.cellId = cellId;
 140       cellData.subId = subId;
 141       cellData.depth = t;
 142       for (i = 0; i < 3; i++) {
 143         cellData.p1[i] = cp1[i];
 144         cellData.p2[i] = cp2[i];
 145       }
 146       /*
 147       cout << "---> Include cell " << cellId << ", depth = " << t   << endl;
 148       cout << "     P1 = (" << cp1[0] << ", " << cp1[1] << ", " << cp1[2] << ")"  << endl;
 149       cout << "     P2 = (" << cp2[0] << ", " << cp2[1] << ", " << cp2[2] << ")"  << endl;
 150       */
 151       IncludeCell(input, cellData, dataList);
 152       if (t < tMin) {
 153         tMin = t;
 154         for (i = 0; i < 3; i++)
 155           xMin[i] = x[i];
 156       }
 157     }
 158   }
 159
 160   if (!dataList.empty()) {
 161     // compare against other actors
 162     IncludeActor(prop3D, input, dataList);
 163     if (tMin < this->GlobalTMin) {
 164       this->MarkPicked(path, prop3D, m, tMin, xMin);
 165     }
 166   }
 167
 168   return tMin;
 169 }
 170
 171 /*!*/
 172 void VTKViewer_CellRectPicker::IncludeCell(vtkDataSet* input,
 173                                            VTKViewer_CellData cellData,
 174                                            VTKViewer_CellDataSet& dataList)
 175 {
 176   vtkGenericCell* cell1 = vtkGenericCell::New();
 177   vtkGenericCell* cell2 = vtkGenericCell::New();
 178   vtkCell* sub1 = 0;
 179   vtkCell* sub2 = 0;
 180
 181   input->GetCell(cellData.cellId, cell1);
 182   if (cell1->GetNumberOfFaces() > 0)
 183     sub1 = cell1->GetFace(cellData.subId);
 184
 185   int i, result;
 186   float p1[3], p2[3], dir[3];
 187   char add = 1;
 188
 189   VTKViewer_CellData curData;
 190   VTKViewer_CellDataSet::iterator it;
 191   for (it = dataList.begin(); it != dataList.end();) {
 192     curData = *it;
 193     for (i = 0; i < 3; i++) {
 194       p1[i] = (cellData.p1[i]+curData.p1[i])/2;
 195       p2[i] = (cellData.p2[i]+curData.p2[i])/2;
 196       dir[i] = p2[i] - p1[i];
 197     }
 198
 199     input->GetCell(curData.cellId, cell2);
 200     sub2 = 0;
 201     if (cell2->GetNumberOfFaces() > 0)
 202       sub2 = cell2->GetFace(curData.subId);
 203
 204     if (sub1) {
 205       if (sub2)
 206         result = IntersectCells(sub1, cellData.depth,
 207         sub2, curData.depth, dir);
 208       else
 209         result = IntersectCells(sub1, cellData.depth,
 210         cell2, curData.depth, dir);
 211     }
 212     else {
 213       if (sub2)
 214         result = IntersectCells(cell1, cellData.depth,
 215         sub2, curData.depth, dir);
 216       else
 217         result = IntersectCells(cell1, cellData.depth,
 218         cell2, curData.depth, dir);
 219     }
 220     /*
 221     cout << "     Direction = (" << dir[0] << ", " << dir[1] << ", " << dir[2] << ")";
 222     cout << "     depth = " << cellData.depth << ", " << curData.depth << ", " << result << endl;
 223     */
 224     if (result > 0) {
 225       if (result == 1) {
 226         it = dataList.erase(it);
 227         continue;
 228       }
 229       else {
 230         add = 0;
 231         break;
 232       }
 233     }
 234     ++it;
 235   }
 236   if (add) {
 237     //cout << "     add "   << endl;
 238     dataList.push_back(cellData);
 239   }
 240
 241   cell1->Delete();
 242   cell2->Delete();
 243 }
 244
 245 /*!*/
 246 void VTKViewer_CellRectPicker::IncludeActor(vtkProp3D* prop,
 247                                             vtkDataSet* input,
 248                                             VTKViewer_CellDataSet& dataList)
 249 {
 250   vtkGenericCell* cell1 = vtkGenericCell::New();
 251   vtkGenericCell* cell2 = vtkGenericCell::New();
 252
 253   int i, result;
 254   float p1[3], p2[3], dir[3];
 255   char removed;
 256
 257   VTKViewer_CellDataSet actorData;
 258   VTKViewer_CellData curData, cellData;
 259   VTKViewer_CellDataSet::iterator it1, it2;
 260   VTKViewer_ActorDataMap::iterator it;
 261
 262   for (it1 = dataList.begin(); it1 != dataList.end();) {
 263     cellData = *it1;
 264     input->GetCell(cellData.cellId, cell1);
 265     removed = 0;
 266
 267     for (it = ActorData.begin(); it != ActorData.end(); ++it) {
 268       //vtkActor* actor = (*it).first;
 269       actorData = (*it).second;
 270
 271       for (it2 = actorData.begin(); it2 != actorData.end();) {
 272         curData = *it2;
 273         for (i = 0; i < 3; i++) {
 274           p1[i] = (cellData.p1[i]+curData.p1[i])/2;
 275           p2[i] = (cellData.p2[i]+curData.p2[i])/2;
 276           dir[i] = p2[i] - p1[i];
 277         }
 278
 279         input->GetCell(curData.cellId, cell2);
 280         result = IntersectCells(cell1, cellData.depth,
 281           cell2, curData.depth, dir);
 282         if (result > 0) {
 283           if (result == 1) {
 284             it2 = actorData.erase(it2);
 285             continue;
 286           }
 287           else {
 288             removed = 1;
 289             it1 = dataList.erase(it1);
 290             break;
 291           }
 292           break;
 293         }
 294         ++it2;
 295       }
 296       if (removed) break;
 297     } // by actors
 298     if (!removed) ++it1;
 299   }
 300
 301   if (!dataList.empty()) {
 302     vtkActor* actor;
 303     if ((actor = vtkActor::SafeDownCast(prop)) != NULL)
 304       ActorData[actor] = dataList;
 305   }
 306
 307   cell1->Delete();
 308   cell2->Delete();
 309 }
 310
 311 /*! The method returns a non-zero value, if the cells intersect each other
 312  * in the direction dir[3].
 313  * \param c1 - first vtkCell pointer
 314  * \param d1 -
 315  * \param c2 - second vtkCell pointer
 316  * \param d2 -
 317  * \param dir[3] - direction
 318  * \retval integer - Returned value is 1 if the first cell is top of
 319  * the second one according to direction and 2 if the second cell is top.
 320  */
 321 int VTKViewer_CellRectPicker::IntersectCells(vtkCell* c1, float d1, vtkCell* c2, float d2, float dir[3])
 322 {
 323   int i, k;
 324   float *orig, cp[3], bounds[6];
 325   vtkCell* cell = 0;
 326   vtkCell* tcell = 0;
 327
 328   // define the central point and cell for projection
 329   if (fabs(d1-d2) < 1.e-7) return 0;
 330
 331   if (d1 < d2) {
 332     orig = c1->Points->GetPoint(0);
 333     cell = c1;
 334     tcell = c2;
 335     c2->GetBounds(bounds);
 336     for (i = 0; i < 3; i++)
 337       cp[i] = (bounds[2*i] + bounds[2*i+1])/2;
 338   }
 339   else if (d1 > d2) {
 340     orig = c2->Points->GetPoint(0);
 341     cell = c2;
 342     tcell = c1;
 343     c1->GetBounds(bounds);
 344     for (i = 0; i < 3; i++)
 345       cp[i] = (bounds[2*i] + bounds[2*i+1])/2;
 346   }
 347
 348   float proj[3];
 349   vtkPlane::GeneralizedProjectPoint(cp, orig, dir, proj);
 350
 351   float *p, pp[3];
 352   vtkPoints* aPoints = vtkPoints::New();
 353   for (k = 0; k < cell->Points->GetNumberOfPoints(); k++) {
 354     p = cell->Points->GetPoint(k);
 355     vtkPlane::GeneralizedProjectPoint(p, orig, dir, pp);
 356     aPoints->InsertNextPoint(pp[0], pp[1], pp[2]);
 357   }
 358
 359   int result = 0;
 360   if (PointInside(proj, aPoints)) {
 361     result = (d1<d2) ? 1 : 2;
 362   }
 363
 364   if (result == 0) {
 365     // test for boundary cells
 366     for (k = 0; k < tcell->Points->GetNumberOfPoints(); k++) {
 367       p = tcell->Points->GetPoint(k);
 368       vtkPlane::GeneralizedProjectPoint(p, orig, dir, pp);
 369       if (PointInside(pp, aPoints)) {
 370         result = 1; break;
 371       }
 372     }
 373     if (result && d2<d1) result = 2;
 374   }
 375
 376   aPoints->Delete();
 377   return result;
 378 }
 379
 380 /*! Check point position.
 381  * \param point[3] - point
 382  * \param list - point data set
 383  * \retval 0 or 1 - The method returns a non-zero value, if the point is inlide point data set
 384  */
 385 char VTKViewer_CellRectPicker::PointInside(float point[3], vtkPoints* list)
 386 {
 387   int i, k, l, n;
 388   float x[3], *a1, *a2, v[3], p[3], q[3];
 389
 390   n = list->GetNumberOfPoints();
 391   if (n < 3) return 0;
 392
 393   for (k = 0; k < n; k++) {
 394     l = k+1; if (l == n) l = 0;
 395     a1 = list->GetPoint(k);
 396     a2 = list->GetPoint(l);
 397
 398     // compute vectors
 399     for (i = 0; i < 3; i++) {
 400       x[i] = point[i] - a1[i];
 401       v[i] = a2[i]-a1[i];
 402     }
 403
 404     // check if point coincides with vertex
 405     if (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] < 1.e-7) return 0;
 406
 407     // compute normal
 408     vtkMath::Cross(x, v, p);
 409
 410     // check if point lays on edge
 411     if (p[0]*p[0] + p[1]*p[1] + p[2]*p[2] < 1.e-7) return 0;
 412
 413     // check if normals are in the same direction
 414     if (k > 0) {
 415       if (vtkMath::Dot(p, q) < 0.0) return 0;
 416     }
 417
 418     // remember normal direction
 419     for (i = 0; i < 3; i++) {
 420       q[i] = p[i];
 421     }
 422   }
 423
 424   return 1;
 425 }