vtkCellLocator::FindcellsAlongLine does not find all the cells

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

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 |
+----+----+----+----+
456 |  7 |
+----+----+----+----+
01 |  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

_______________________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://vtk.org/mailman/listinfo/vtkusers