[vtkusers] Minor problem with vtkCellLocator::IntersectWithLine

[Paul A Hsieh]

The method vtkCellLocator::IntersectWithLine appears to not work correctly
under the following situation:
1. When the data set is a vtkPolyData that is entirely flat and lies in the
plane z = constant,
AND 2. The sampling line is parallel to the z axis.
This problem is illustrated by the attached C++ sample program (end of this
email). The poly data consist of a single cell, which is a square that lies
on the plane z = 0 and centered on the z axis. The line runs along the z
axis from z = -10 to z = 10. This line should intersect the cell. However,
vtkCellLocator::IntersectWithLine returns 0. If the cell is slightly
tilted, then vtkCellLocator::IntersectWithLine correctly returns 1. This
occurred with both vtk 4.0 and the latest nightly release (Nov 23, 2002). I
am running on Windows 2000.

Looking at the source code (snippet below) it would appear that for the
above case, tmax would end up being zero at the end of the first "for"
loop. This would cause a division by zero in the second "for" loop.
Paul Hsieh



int vtkCellLocator::IntersectWithLine(float a0[3], float a1[3], float tol,
                                      float& t, float x[3], float
pcoords[3],
                                      int &subId, vtkIdType &cellId,
                                      vtkGenericCell *cell)
{

  /* ... skipped ... */

  // get the bounds
  bounds = this->DataSet->GetBounds();

  // convert the line into i,j,k coordinates
  tMax = 0.0;
  for (i=0; i < 3; i++)
    {
    direction1[i] = a1[i] - a0[i];
    if ( (bounds[2*i+1] - bounds[2*i]) != 0.0)
      {
      origin[i] = (a0[i] - bounds[2*i])/(bounds[2*i+1] - bounds[2*i]);
      direction2[i] = (a1[i] - a0[i])/(bounds[2*i+1] - bounds[2*i]);
      }
    else
      {
      origin[i] = 0.0;
      direction2[i] = 0.0;
      }
    bounds2[2*i]   = 0.0;
    bounds2[2*i+1] = 1.0;
    tMax += direction2[i]*direction2[i];
    }

  tMax = sqrt(tMax);
  stopDist = tMax*this->NumberOfDivisions;
  for (i = 0; i < 3; i++)
    {
    direction3[i] = direction2[i]/tMax;
    }


============Sample code to illustrate the problem ===========


#include "vtkPolyData.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkCellLocator.h"

void main()
{
  int i;
  static float x[4][3]={{-1,-1,0}, {1,-1,0}, {1,1,0}, {-1,1,0}};
  static vtkIdType pts[4]={0,1,2,3};

  vtkPolyData *pd = vtkPolyData::New();
  vtkPoints *points = vtkPoints::New();
  vtkCellArray *polys = vtkCellArray::New();

  for (i=0; i<4; i++) points->InsertPoint(i,x[i]);
  for (i=0; i<1; i++) polys->InsertNextCell(4,pts);

  pd->SetPoints(points);
  points->Delete();
  pd->SetPolys(polys);
  polys->Delete();

  vtkCellLocator *cellLocator = vtkCellLocator::New();
  cellLocator->SetDataSet(pd);
  cellLocator->BuildLocator();

  float p1[] = {1, 0, 10};
  float p2[] = {1, 0, -10};
  float t, ptline[3], pcoords[3];
  int subId;
  int result = cellLocator->IntersectWithLine(p1, p2, 0.001, t, ptline,
pcoords, subId);

  // "result" should be 1 but the above method returns zero.
  // However, if the cell is tilted by displacing one of the points
  // slightly higher, then the method return 1.
  cout << result << endl;

  cellLocator->Delete();
  pd->Delete();
}