[vtkusers] vtkCellLocator::FindcellsAlongLine does not find all the cells

[Alexander Pletzer]

Hi VTK'ians,

I found vtkCellLocator to be a very fast way to locate  cells but I now
stumbled on an issue where the cell locator does not locate all the cells
along a line. In the code below, an unstructured 2x4x4 grid and
corresponding vtkCellLocator objects are created. The grid is a collection
of cells. A line is created that runs along the the bottom of the grid,
crossing all the latitudes and longitudes. Shown below are the bottom cells:

+----+----+----+----+
| 12 | 13 | *14* | *15* |
+----+----+----+----+
|  8 |  *9* | *10* | *11* |
+----+----+----+----+
|  *4* |  *5* |  *6* |  7 |
+----+----+----+----+
|  *0* |  *1* |  2 |  3 |
+----+----+----+----+

I would expect the number of cells found by FindCellsAlongLine to include
0, 5, 10 and 15 (red cells) as these cells are crossed by the line. Cells
1, 4, 6, 9, 11 and 14 touch the line so should also be found (blue cells).
The total number of cells is 10. When running the code (VTK 8.1.1 on mac) I
get:

...

1: line intersects with cell 0

1: line intersects with cell 1

1: line intersects with cell 4

1: line intersects with cell 5

1: line intersects with cell 6

1: line intersects with cell 9

1: line intersects with cell 10
1: line intersects with cell 14

 We see that cells 11 and 15 are missing! Other resolutions return the
correct list. I played with the tolerance in FindCellsAlongLine but this
seems to have no effect.

Any idea what I'm doing wrong? Thanks in advance.

--Alex

Test code:

#undef NDEBUG // turn on asserts
#include <vtkUnstructuredGrid.h>
#include <vtkCellLocator.h>
#include <vtkPoints.h>
#include <vtkIdList.h>
#include <vtkGenericCell.h>
#include <iostream>
#include <cmath>

void testLatLon(size_t nElv, size_t nLat, size_t nLon, const double pa[],
const double pb[]) {

    //
    // generate vertices
    //
    double p[3];

    double dElv = 1.0 / double(nElv);
    double dLat = 180.0 / double(nLat);
    double dLon = 360.0 / double(nLon);

    vtkPoints* points = vtkPoints::New();
    size_t nCells = nElv * nLat * nLon;
    points->SetDataTypeToDouble();
    points->SetNumberOfPoints(8 * nCells);

    vtkIdType index = 0;
    vtkIdType cellId = 0;
    for (size_t k = 0; k < nElv; ++k) {
        double elv0 = dElv*k;
        double elv1 = elv0 + dElv;
        for (size_t j = 0; j < nLat; ++j) {
            double lat0 = -90.0 + j*dLat;
            double lat1 = lat0 + dLat;
            for (size_t i = 0; i < nLon; ++i) {
                double lon0 = -180.0 + i*dLon;
                double lon1 = lon0 + dLon;

                p[0] = elv0; p[1] = lat0; p[2] = lon0;
                points->SetPoint(index, p);
                index++;

                p[0] = elv0; p[1] = lat0; p[2] = lon1;
                points->SetPoint(index, p);
                index++;

                p[0] = elv0; p[1] = lat1; p[2] = lon1;
                points->SetPoint(index, p);
                index++;

                p[0] = elv0; p[1] = lat1; p[2] = lon0;
                points->SetPoint(index, p);
                index++;

                p[0] = elv1; p[1] = lat0; p[2] = lon0;
                points->SetPoint(index, p);
                index++;

                p[0] = elv1; p[1] = lat0; p[2] = lon1;
                points->SetPoint(index, p);
                index++;

                p[0] = elv1; p[1] = lat1; p[2] = lon1;
                points->SetPoint(index, p);
                index++;

                p[0] = elv1; p[1] = lat1; p[2] = lon0;
                points->SetPoint(index, p);
                index++;

                std::cout << "cell " << cellId << ": elv = " << elv0 << "
-> " << elv1 <<
                             " lat = " << lat0 << " -> " << lat1 <<
                             " lon = " << lon0 << " -> " << lon1 << '\n';

                cellId++;
            }
        }
    }

    //
    // create grid
    //
    vtkUnstructuredGrid* grid = vtkUnstructuredGrid::New();
    grid->SetPoints(points);
    grid->Allocate(1, 1);
    vtkIdList* ptIds = vtkIdList::New();
    ptIds->SetNumberOfIds(8);
    for (size_t iCell = 0; iCell < nCells; ++iCell) {
        for (vtkIdType iVert = 0; iVert < 8; ++iVert) {
            ptIds->SetId(iVert, 8*iCell + iVert);
        }
        grid->InsertNextCell(VTK_HEXAHEDRON, ptIds);
        std::cout << "inserting cell " << iCell;
        for (vtkIdType iVert = 0; iVert < 8; ++iVert) {
          vtkIdType vi = ptIds->GetId(iVert);
          double vrtx[3];
          points->GetPoint(vi, vrtx);
          std::cout << "\t" << iVert << " vert index: " << vi <<
          " vert coords: " << vrtx[0] << ',' << vrtx[1] << ',' << vrtx[2]
<< '\n';
        }
    }

    vtkCellLocator* loc = vtkCellLocator::New();
    loc->SetDataSet(grid);
    loc->BuildLocator();

    // find all the cells allong the line
    vtkIdList* cellIds = vtkIdList:: New();
    double tol = 0.01;
    loc->FindCellsAlongLine((double*) pa, (double*) pb, tol, cellIds);
    for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i) {
      std::cout << "line intersects with cell " << cellIds->GetId(i) <<
'\n';
    }

    // clean up
    cellIds->Delete();
    loc->Delete();
    ptIds->Delete();
    grid->Delete();
    points->Delete();
}


int main(int argc, char** argv) {

    double pa[] = {0., -90., -180.};
    double pb[] = {0., 90., 180.};
    testLatLon(2, 4, 4, pa, pb);

    return 0;
}






-- 
Alexander Pletzer
HPC Scientific Programmer for New Zealand eScience Infrastructure (NeSI),
Wellington NZ
mobile: +64 21 085 79661
http://www. <http://www.niwa.co.nz>nesi.org.nz
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://vtk.org/pipermail/vtkusers/attachments/20190108/bba5b4ee/attachment.html>