How to get texture coordinates from VTK IntersectWithLine?

Question

I've loaded a texture mapped OBJ via vtkOBJReader and loaded it into a vtkModifiedBSPTree:

auto readerOther(vtkSmartPointer::New());
auto rawOtherPath(modelPathOther.toLatin1());
readerOther->SetFileName(rawOtherPath.data());
readerOther->Update();

auto meshDataOther(readerOther->GetOutput());

auto bspTreeOther(vtkSmartPointer::New());
bspTreeOther->SetDataSet(meshDataOther);
bspTreeOther->BuildLocator();

I then compute my line segment start and end and feed that into

if (bspTreeOther->IntersectWithLine(p1, p2, tolerance, distanceAlongLine, intersectionCoords, pcoords, subId, cellId, cell))

With all the relevant predefined variables of course.

What I need is the texture's UV coordinates at the point of intersection.

I'm so very new to VTK that I've not yet caught the logic of how its put together yet; the abstraction layers are still losing me while I'm digging through the source.

I've hunted for this answer across SO and the VTK users archives and found vague hints given by those who understood VTK deeply to those who were nearly there themselves, and thus of little help to me thus far.

(Appended 11/9/2018) To clarify, I'm working with non-degenerate triangulated meshes created by a single 3D scanner shot, so quads and other higher polygons are not going to be ever seen by my code. A general solution should account for such things, but that can be accomplished via triangulating the mesh first via a good application of handwavium.

Charlie · Accepted Answer

Code

Note that if one vertex belongs to several polygons and has different texture coordinates, VTK will create duplicates of the vertex. I don't use vtkCleanPolyData, because VTK will merge such "duplicates" and we will lose needed information, as far as I know.

I use vtkCellLocator instead of vtkModifiedBSPTree, because in my case it was faster.

The main file main.cpp. You can find magic numbers in start and end arrays — these are your p1 and p2. I've set these values just for example

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 

int main(int argc, char * argv[])
{
    if (argc < 2)
    {
        std::cerr << "Usage: " << argv[0] << " OBJ_file_name" << std::endl;
        return EXIT_FAILURE;
    }
    auto reader{vtkSmartPointer::New()};
    reader->SetFileName(argv[1]);
    reader->Update();
    // Triangulate the mesh if needed
    auto triangleFilter{vtkSmartPointer::New()};
    triangleFilter->SetInputConnection(reader->GetOutputPort());
    triangleFilter->Update();
    auto mesh{triangleFilter->GetOutput()};
    // Use `auto mesh(reader->GetOutput());` instead if no triangulation needed

    // Build a locator to find intersections
    auto locator{vtkSmartPointer::New()};
    locator->SetDataSet(mesh);
    locator->BuildLocator();

    // Initialize variables needed for intersection calculation
    double start[3]{-1, 0, 0.5};
    double end[3]{   1, 0, 0.5};
    double tolerance{1E-6};
    double relativeDistanceAlongLine;
    double intersectionCoordinates[3];
    double parametricCoordinates[3];
    int subId;
    vtkIdType cellId;
    auto cell{vtkSmartPointer::New()};
    // Find intersection
    int intersected = locator->IntersectWithLine(
        start,
        end,
        tolerance,
        relativeDistanceAlongLine,
        intersectionCoordinates,
        parametricCoordinates,
        subId,
        cellId,
        cell.Get()
    );

    // Get points of intersection cell
    auto pointsIds{vtkSmartPointer::New()};
    mesh->GetCellPoints(cellId, pointsIds);
    // Store coordinates and texture coordinates of vertices of the cell
    double meshTrianglePoints[3][3];
    double textureTrianglePoints[3][2];
    auto textureCoordinates{mesh->GetPointData()->GetTCoords()};
    for (unsigned pointNumber = 0; pointNumber < cell->GetNumberOfPoints(); ++pointNumber)
    {
        mesh->GetPoint(pointsIds->GetId(pointNumber), meshTrianglePoints[pointNumber]);
        textureCoordinates->GetTuple(pointsIds->GetId(pointNumber), textureTrianglePoints[pointNumber]);
    }

    // Normalize the coordinates
    double movedMeshTrianglePoints[3][3];
    for (unsigned i = 0; i < 3; ++i)
    {
        movedMeshTrianglePoints[0][i] = 0;
        movedMeshTrianglePoints[1][i] =
            meshTrianglePoints[1][i] -
            meshTrianglePoints[0][i];
        movedMeshTrianglePoints[2][i] =
            meshTrianglePoints[2][i] -
            meshTrianglePoints[0][i];
    }
    // Normalize the texture coordinates
    double movedTextureTrianglePoints[3][2];
    for (unsigned i = 0; i < 2; ++i)
    {
        movedTextureTrianglePoints[0][i] = 0;
        movedTextureTrianglePoints[1][i] =
            textureTrianglePoints[1][i] -
            textureTrianglePoints[0][i];
        movedTextureTrianglePoints[2][i] =
            textureTrianglePoints[2][i] -
            textureTrianglePoints[0][i];
    }

    // Calculate SVD of a matrix consisting of normalized vertices
    double U[3][3];
    double w[3];
    double VT[3][3];
    vtkMath::SingularValueDecomposition3x3(movedMeshTrianglePoints, U, w, VT);
    // Calculate pseudo inverse of a matrix consisting of normalized vertices
    double pseudoInverse[3][3]{0};
    for (unsigned i = 0; i < 3; ++i)
    {
        for (unsigned j = 0; j < 3; ++j)
        {
            for (unsigned k = 0; k < 3; ++k)
            {
                if (w[k] != 0)
                {
                    pseudoInverse[i][j] += VT[k][i] * U[j][k] / w[k];
                }
            }
        }
    }

    // Calculate interpolation matrix
    double interpolationMatrix[3][2]{0};
    for (unsigned i = 0; i < 3; ++i)
    {
        for (unsigned j = 0; j < 2; ++j)
        {
            for (unsigned k = 0; k < 3; ++k)
            {
                interpolationMatrix[i][j] += pseudoInverse[i][k] * movedTextureTrianglePoints[k][j];
            }
        }
    }

    // Calculate interpolated texture coordinates of the intersection point
    double interpolatedTexturePoint[2]{textureTrianglePoints[0][0], textureTrianglePoints[0][1]};
    for (unsigned i = 0; i < 2; ++i)
    {
        for (unsigned j = 0; j < 3; ++j)
        {
            interpolatedTexturePoint[i] += (intersectionCoordinates[j] - meshTrianglePoints[0][j]) * interpolationMatrix[j][i];
        }
    }

    // Print the result
    std::cout << "Interpolated texture coordinates";
    for (unsigned i = 0; i < 2; ++i)
    {
        std::cout << " " << interpolatedTexturePoint[i];
    }
    std::cout << std::endl;

    return EXIT_SUCCESS;
}

CMake project file CMakeLists.txt

cmake_minimum_required(VERSION 3.1)

PROJECT(IntersectInterpolate)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

add_executable(IntersectInterpolate MACOSX_BUNDLE main.cpp)

if(VTK_LIBRARIES)
    target_link_libraries(IntersectInterpolate ${VTK_LIBRARIES})
else()
    target_link_libraries(IntersectInterpolate vtkHybrid vtkWidgets)
endif()

Math

What we need

Suppose you have a mesh consisting of triangles and your vertices have texture coordinates.

Given vertices of a triangle A, B and C, corresponding texture coordinates A', B' and C', you want to find a mapping (to interpolate) from another inner and boundary points of the triangle to the texture. Let's make some rational assumptions:

Points A, B, C should correspond to their texture coordinates A', B', C';
Each point X on the border, say AB, should correspond to the points of A'B' line in following way: |AX| / |AB| = |A'X'| / |A'B'| — half way on the original triangle should be a half way on the texture map;
- Centroid of the triangle (A + B + C) / 3 should correspond to centroid of the texture triangle (A' + B' + C') / 3.

Equations to solve

Looks like we want to have affine mapping: coordinates of vertices of the original triangle should be multiplied by some coefficients and be added to some constants. Let's construct the system of equations

Ax * Mxx + Ay * Myx + Az * Mzx + M0x = A'x
Ax * Mxy + Ay * Myy + Az * Mzy + M0y = A'y
Ax * Mxz + Ay * Myz + Az * Mzz + M0z = 0

and the same for B and C. You can see that we have 9 equations and 12 unknowns. Though, equations containing Miz (for i in {x, y, z}) have the solution 0 and don't play any role in further computations, so we can just set them equal to 0. Thus, we have system with 6 equations and 8 unknowns

Ax * Mxx + Ay * Myx + Az * Mzx + M0x = A'x
Ax * Mxy + Ay * Myy + Az * Mzy + M0y = A'y

Let's write entire system in matrix view

--          --   --       --   --       --
| 1 Ax Ay Az |   | M0x M0y |   | A'x A'y |
| 1 Bx By Bz | x | Mxx Mxy | = | B'x B'y |
| 1 Cx Cy Cz |   | Myx Myy |   | C'x C'y |
--          --   | Mzx Mzy |   --       --
                 --       --

I subtract coordinates of A vertex from B and C and texture coordinates A' from B' and C', and now we have the triangle with the first vertex located in start of coordinate system as well as corresponding texture coordinates. This means that now triangles are not translated (moved) one relative to another and we don't need M0 part of interpolation matrix

--        --   --       --   --       --
| Bx By Bz |   | Mxx Mxy |   | B'x B'y |
| Cx Cy Cz | x | Myx Myy | = | C'x C'y |
--        --   | Mzx Mzy |   --       --
               --       --

Solution

Let's call the first matrix P, the second M and the last one T

P M = T

The matrix P is not square. If we add zero row to it, the matrix becomes singular. So, we have to calculate pseudo-inverse of it in order to solve the equation. There's no function for calculating pseudo-inverse matrix in VTK. We go to Moore–Penrose inverse article on Wikipedia and see that it can be calculated using SVD. VTKMath::SingularValueDecomposition3x3 function allows us to do it. The function gives us U, S and VT matrices. I'll write pseudo-inverse of matrix P as P", transposition of U as UT and transposition of VT as V. Pseudo-inverse of diagonal matrix S is a matrix with 1 / Sii elements where Sii is not a zero and 0 for zero elements

P = U S VT
P" = V S" UT
M = P" T

Usage

To apply interpolation matrix, we need to not forget that we need to translate input and output vectors. A' is a 2D vector of texture coordinates of the first vertex in the triangle, A is a 3D vector of coordinates of the vertex, M is the found interpolation matrix, p is a 3D intersection point we want to get texture coordinates for, t' is the resulting 2D vector with interpolated texture coordinates

t' = A' + (p - A) M