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