"""
\file test_vtkImplicitDataSet_tol2_problem.py

\date 2011-07-04

\author Peter Wood (peter@irba.co.nz)

\brief VTK version 5.6.1, 

Tolerance problem in 'vtkImplicitDataSet.cxx' in method
'EvaluateFunction' line 81

// Find the cell that contains xyz and get it
cell = this->DataSet->FindAndGetCell(x,NULL,-1,0.0,subId,pcoords,this->Weights);
                                               ^^^
Note 4th parameter '0.0'.  This is 'tol2' in API, which is used in
'vtkPointSet.cxx' in method 'FindCellWalk'

    double dist2;
    if (   (cell->EvaluatePosition(x, closestPoint, subId,
                                   pcoords, dist2, weights) == 1)
        && (dist2 <= tol2) )
      {
      return cellId;
      }

The problem is that 'tol2' is 0.0, but often 'dist2' may be calculated
to be approx 1e-032 to 1e-029, a very small number, but still greater
than 0.0! (NB: I adapted vtkPointSet.cxx to report 'dist2')

To fix this problem, when calling 'FindAndGetCell' in
'vtkImplicitDataSet' replace '0.0' with an appropriate small positive
tolerance value, e.g. '1e-16'.  This could be user-configurable?

"""

import numpy as np
import vtk

NUM_INTERP_POINTS = 100

num_pts = 4

pts = vtk.vtkPoints()
pts.SetNumberOfPoints(num_pts)
pts.InsertPoint(0, (0.0, 5.0, 0.0))
pts.InsertPoint(1, (10.0, 5.0, 0.0))
pts.InsertPoint(2, (20.0, 5.0, 0.0))
pts.InsertPoint(3, (30.0, 5.0, 0.0))

cell_array = vtk.vtkCellArray()

# three lines between the points
for i in range(3):
    line = vtk.vtkLine()
    pt_ids = line.GetPointIds()
    pt_ids.SetId(0, i)
    pt_ids.SetId(1, i+1)
    cell_array.InsertNextCell(line)
#
# create data array for points data
z_array = vtk.vtkDoubleArray()
z_array.SetName("z_height")
z_array.SetNumberOfComponents(1)
z_array.Allocate(num_pts, 2)
z_array.InsertValue(0, 0.0)
z_array.InsertValue(1, 0.0)
z_array.InsertValue(2, 5.0)
z_array.InsertValue(3, 5.0)

polydata = vtk.vtkPolyData()
polydata.SetPoints(pts)
polydata.SetLines(cell_array)
pt_data = polydata.GetPointData()
pt_data.SetScalars(z_array)

_out_val = -99.0
implicit_func = vtk.vtkImplicitDataSet()
implicit_func.SetDataSet(polydata)
implicit_func.SetOutValue(_out_val)
#
pt_start_np = np.array((0.0, 5.0, 0.0))
ilen = 30.0
dir_v = np.array((1.0, 0.0, 0.0))
eps = 1e-6

# 1D interpolation on scalars 'z_array' along polyline.
for i in range(NUM_INTERP_POINTS + 1):
    i_pt = pt_start_np + ((float(i) * ilen) / float(NUM_INTERP_POINTS)) * dir_v

    _z = implicit_func.EvaluateFunction(tuple(i_pt))
    print tuple(i_pt), _z,
    if abs(_out_val - _z) < eps:
        # _z is OutVale
        print '<***** problem '
    else:
        print ''