[vtkusers] vtkCellLocator precision issue

[xabi riobe]

Hi,

I have found a numerical precision problem with the
method vtkCellLocator::IntersectWithLine.

in line 287 of revision 1.89:
  if (cell->IntersectWithLine(a0, a1, tol, t, x, pcoords, subId) )
    {
    if ( ! this->IsInOctantBounds(x) )

sometimes the resulting intersection point returned by the
cell->IntersectWithLine is a (0.000000000000004) little out of the bounds
and the call to IsInOctantBounds discard the point :(

following is an example that illustrates that.
a vtkCellLocator is filled with a cube and a line going from dPoint to
dPoint2 intersect the box.
in Z axis, the bounds of the box are (1.0; -20)

if dPoint2[2] = -77.76333863288239, the resulting intersection point is at
-20 and all is correct
if dPoint2[2] = -77.76333863288238, the resulting intersection point is
at -20.000000000000004, so out of bounds.


So will it be a good thing to add a tolerance to the
method vtkCellLocator::IsInOctantBounds(double x[3]) ?



#include <vtk/vtkActor.h>
#include <vtk/vtkCamera.h>
#include <vtk/vtkCellLocator.h>
#include <vtk/vtkCubeSource.h>
#include <vtk/vtkGenericCell.h>
#include <vtk/vtkInteractorStyleTrackballCamera.h>
#include <vtk/vtkLineSource.h>
#include <vtk/vtkPolyData.h>
#include <vtk/vtkPolyDataMapper.h>
#include <vtk/vtkProperty.h>
#include <vtk/vtkRenderer.h>
#include <vtk/vtkRenderWindow.h>
#include <vtk/vtkRenderWindowInteractor.h>
#include <vtk/vtkSmartPointer.h>
#include <vtk/vtkSphereSource.h>
#include <vtk/vtkSTLReader.h>
#include <vtk/vtkTriangleFilter.h>

#define MY_SP(class, variable)\
   vtkSmartPointer<class> variable = vtkSmartPointer<class>::New();

double dPoint[3] = {
   19.277858734130859,
   23.247093200683594,
   -0.45018148422241211
};

double dPoint2[3] = {
   19.277858734130859,
   23.247093200683594,
   // error
   -77.76333863288238
   // ok
   //-77.76333863288239
};

//---------------------------------------------------------------------------------------------
void Test(vtkPolyData *data, double sourcePnt[3], double destinPnt[3])
{
   MY_SP(vtkCellLocator, loc);
   loc->SetTolerance(0.0);
   loc->FreeSearchStructure();
   loc->CacheCellBoundsOn();
   loc->AutomaticOn();
   loc->SetDataSet(data);
   loc->BuildLocator();
   loc->Update();

   vtkIdType cellId;
   int subId;
   double param_t, intersect[3], paraCoord[3];
   vtkGenericCell *cell = vtkGenericCell::New();

   if(! loc->IntersectWithLine(sourcePnt, destinPnt, 0.0010, param_t,
      intersect, paraCoord, subId, cellId, cell) )
   {
      cout << "error" << endl;
   }

   cout << "intersection point: " << intersect[0] << " " << intersect[1] <<
" " << intersect[2] << endl;

   cell->Delete();
}

//---------------------------------------------------------------------------------------------
int main(int , char* [])
{
   MY_SP(vtkRenderer, ren1);
   ren1->SetBackground(0.2, 0.2, 0.2);
   MY_SP(vtkRenderWindow, renWin);
   renWin->SetSize( 800, 800 );
   renWin->AddRenderer(ren1);
   MY_SP(vtkRenderWindowInteractor, iren);
   iren->SetRenderWindow(renWin);

   MY_SP(vtkInteractorStyleTrackballCamera, vtk_style);
   iren->SetInteractorStyle(vtk_style);

   /////
   MY_SP(vtkCubeSource, cube);
   cube->SetBounds(-43.5, 43.5, -35.0, 35.0, -20.0, 1.0);

   MY_SP(vtkTriangleFilter, tri);
   tri->SetInputConnection(cube->GetOutputPort());
   tri->Update();

   MY_SP(vtkPolyDataMapper, Mapper);
   Mapper->SetInputConnection(tri->GetOutputPort());
   Mapper->ImmediateModeRenderingOn();

   MY_SP(vtkActor, Actor);
   Actor->SetMapper(Mapper);
   Actor->GetProperty()->SetRepresentationToWireframe();

   ren1->AddActor(Actor);

   //////
   MY_SP(vtkSphereSource, sph);
   sph->SetCenter(dPoint);
   sph->SetRadius(1);
   MY_SP(vtkPolyDataMapper, MapperSph);
   MapperSph->SetInputConnection(sph->GetOutputPort());
   MapperSph->ImmediateModeRenderingOn();

   MY_SP(vtkActor, ActorSph);
   ActorSph->SetMapper(MapperSph);
   ActorSph->GetProperty()->SetColor(1.0, 0.0, 0.0);
   ren1->AddActor(ActorSph);

   MY_SP(vtkSphereSource, sph2);
   sph2->SetCenter(dPoint2);
   sph2->SetRadius(1);
   MY_SP(vtkPolyDataMapper, MapperSph2);
   MapperSph2->SetInputConnection(sph2->GetOutputPort());
   MapperSph2->ImmediateModeRenderingOn();

   MY_SP(vtkActor, ActorSph2);
   ActorSph2->SetMapper(MapperSph2);
   ActorSph2->GetProperty()->SetColor(0.0, 0.0, 1.0);
   ren1->AddActor(ActorSph2);

   MY_SP(vtkLineSource, line);
   line->SetPoint1(dPoint);
   line->SetPoint2(dPoint2);
   MY_SP(vtkPolyDataMapper, Mapperline);
   Mapperline->SetInputConnection(line->GetOutputPort());
   Mapperline->ImmediateModeRenderingOn();

   MY_SP(vtkActor, Actorline);
   Actorline->SetMapper(Mapperline);
   Actorline->GetProperty()->SetColor(0.0, 1.0, 0.0);
   ren1->AddActor(Actorline);

   Test(cube->GetOutput(), dPoint, dPoint2);

   ren1->ResetCamera();
   renWin->Render();
   iren->Start();

   ren1 = 0;

   return 0;
}
//---------------------------------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20100204/464ff851/attachment.htm>