results of Curved PN-Triangles algorithm has visible edges

Question

I implemented the curved PN triangles algorithm described in this paper to smooth geometry in real time. After the implementation was finished (in glsl) I realized, that this algorithm smooth the geometry but there are still edges visible. To better debug I converted the program to C++ and inserted a unit cube. Important: the unit cubes faces share their vertices... so we have 8 vertices with 8 normal vectors (no flat shading).

The used obj file (cube.obj):

v -1 -1 -1
v 1 -1 -1
v -1 1 -1
v 1 1 -1
v -1 -1 1
v 1 -1 1
v -1 1 1
v 1 1 1
vn -1 -1 -1
vn 1 -1 -1
vn -1 1 -1
vn 1 1 -1
vn -1 -1 1
vn 1 -1 1
vn -1 1 1
vn 1 1 1
f 1//1 2//2 3//3
f 3//3 2//2 4//4
f 2//2 6//6 4//4
f 4//4 6//6 8//8
f 3//3 4//4 7//7
f 7//7 4//4 8//8
f 1//1 5//5 2//2
f 2//2 5//5 6//6
f 1//1 3//3 5//5
f 5//5 3//3 7//7
f 5//5 7//7 6//6
f 6//6 7//7 8//8

as you can see the normal vectors (vn) are not normalized. This is no problem, because they will be normalized during filling the buffer.

The following code reads the mesh file and subdivides each face with respect to the edgeDivisionCount variable. Here it is set to 100. That means each edge is subdivided into 100 vertices.

Code:

#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <iostream>
#include <thread>

struct Vector3D
{
public:
    Vector3D()
    {
        x = 0;
        y = 0;
        z = 0;
    }

    Vector3D(float x, float y, float z)
    {
        this->x = x;
        this->y = y;
        this->z = z;
    }

    Vector3D operator*(float scalar) const
    {
        return Vector3D(x * scalar, y * scalar, z * scalar);
    }

    Vector3D operator+(const Vector3D& v2) const
    {
        return Vector3D(x + v2.x, y + v2.y, z + v2.z);
    }

    Vector3D operator-(const Vector3D& v2) const
    {
        return Vector3D(x - v2.x, y - v2.y, z - v2.z);
    }

    Vector3D operator-() const
    {
        return Vector3D(-x, -y, -z);
    }

    float dot(const Vector3D& v2) const
    {
        return x*v2.x + y*v2.y + z*v2.z;
    }

    float length() const
    {
        return sqrt(x*x + y*y + z*z);
    }

    Vector3D normalize() const
    {
        float l = length();
        return Vector3D(x / l, y / l, z / l);
    }

    float x;
    float y;
    float z;
};

void saveOBJ(std::string file, std::vector<Vector3D>* vertices, std::vector<std::tuple<int, int, int>>* faces)
{
    std::ofstream myfile;
    myfile.open(file);  
    if (vertices != nullptr)
    {
        for (std::vector<Vector3D>::iterator it = vertices->begin(); it != vertices->end(); ++it)
        {
            myfile << "v " << it->x << " " << it->y << " " << it->z << "\n";
        }
    }
    if (faces != nullptr)
    {
        for (std::vector<std::tuple<int, int, int>>::iterator it = faces->begin(); it != faces->end(); ++it)
        {
            myfile << "f " << std::get<0>(*it) + 1 << " " << std::get<1>(*it) + 1 << " " << std::get<2>(*it) + 1 << "\n";
        }
    }
    myfile.close();
}

Vector3D curvedPNVertex(Vector3D b300, Vector3D b030, Vector3D b003, Vector3D n200, Vector3D n020, Vector3D n002, Vector3D b210, Vector3D b120, Vector3D b021, Vector3D b012, Vector3D b102, Vector3D b201, Vector3D b111, float scalar0, float scalar1, float scalar2)
{
    float u = scalar1;
    float v = scalar2;
    float w = scalar0;
    float u2 = pow(u,2);
    float v2 = pow(v,2);
    float w2 = pow(w,2);
    float u3 = pow(u,3);
    float v3 = pow(v,3);
    float w3 = pow(w,3);

    return b003*v3 + b102*3*w*v2 + b012*3*u*v2 + b201*3*w2*v + b111*6*w*u*v + b021*3*u2*v + b300*w3 + b210*3*w2*u + b120*3*w*u2 + b030*u3;
}

std::vector<Vector3D> curvedPNTriangle(int divisionCount, Vector3D v0, Vector3D v1, Vector3D v2, Vector3D n0, Vector3D n1, Vector3D n2)
{
    std::vector<Vector3D> result;

    Vector3D b300 = v0;
    Vector3D b030 = v1;
    Vector3D b003 = v2;

    Vector3D n200 = n0.normalize();
    Vector3D n020 = n1.normalize();
    Vector3D n002 = n2.normalize();

    Vector3D b210 = (b300 * 2.0f + b030 - (n200 * (b030 - b300).dot(n200))) * (1.0f / 3.0f);
    Vector3D b120 = (b030 * 2.0f + b300 - (n020 * (b300 - b030).dot(n020))) * (1.0f / 3.0f);
    Vector3D b021 = (b030 * 2.0f + b003 - (n020 * (b003 - b030).dot(n020))) * (1.0f / 3.0f);
    Vector3D b012 = (b003 * 2.0f + b030 - (n002 * (b030 - b003).dot(n002))) * (1.0f / 3.0f);
    Vector3D b102 = (b003 * 2.0f + b300 - (n002 * (b300 - b003).dot(n002))) * (1.0f / 3.0f);
    Vector3D b201 = (b300 * 2.0f + b003 - (n200 * (b003 - b300).dot(n200))) * (1.0f / 3.0f);
    Vector3D e = (b210+b120+b201+b102+b012+b021) * (1.0f / 6.0f);
    Vector3D b111 = e + (e - (b300+b030+b003) * (1.0f / 3.0f)) * (1.0f / 2.0f);

    //divide triangle face
    for (unsigned int i = 0; i < divisionCount; ++i)
    {
        for (unsigned int j = 0; j < divisionCount - i; ++j)
        {
            float scalar1 = static_cast<float>(i) / divisionCount;
            float scalar0 = static_cast<float>(j) / divisionCount;
            float scalar2 = 1 - (scalar0 + scalar1);
            Vector3D vertex = curvedPNVertex(b300, b030, b003, n200, n020, n002, b210, b120, b021, b012, b102, b201, b111, scalar0, scalar1, scalar2);
            result.push_back(vertex);
        }
    }

    return result;
}

std::vector<int> getIntOutOfString(std::string text)
{
    std::vector<int> result;

    std::string subtext = "";
    for (unsigned int i = 0; i < text.size(); ++i)
    {
        if ((text.at(i) >= '0' && text.at(i) <= '9'))
        {
            subtext += text.at(i);
        }
        else if (subtext.size() == 0 && text.at(i) == '-')
        {
            subtext += text.at(i);
        }
        else if (subtext.size() > 0)
        {
            //end of a int reached! Save int into list!
            result.push_back(std::stoi(subtext));
            subtext = "";
        }
    }
    if (subtext != "")
    {
        result.push_back(std::stoi(subtext));
    }
    return result;
}


std::vector<double> getDoublesOutOfString(std::string text)
{
    std::vector<double> result;

    std::string subtext = "";
    for (unsigned int i = 0; i < text.size(); ++i)
    {
        if ((text.at(i) >= '0' && text.at(i) <= '9') || text.at(i) == '.')
        {
            subtext += text.at(i);
        }
        else if (subtext.size() == 0 && text.at(i) == '-')
        {
            subtext += text.at(i);
        }
        else if (subtext.size() > 0)
        {
            //end of a double reached! Save double into list!
            result.push_back(std::stod(subtext));
            subtext = "";
        }
    }
    if (subtext != "")
    {
        result.push_back(std::stod(subtext));
    }
    return result;
}

bool readOBJ(std::string file, std::vector<Vector3D>* positions, std::vector<Vector3D>* normals, std::vector<std::tuple<int, int, int, int, int, int>>* faces)
{
    std::string line;
    std::ifstream myfile(file);
    if (myfile.is_open())
    {
        while (std::getline(myfile, line))
        {
            if (line.size() > 1)
            {
                if (line[0] == 'v' && line[1] == ' ')
                {
                    std::vector<double> vertexValues = getDoublesOutOfString(line);
                    if (vertexValues.size() == 3)
                    {
                        Vector3D position = Vector3D(vertexValues[0], vertexValues[1], vertexValues[2]);
                        positions->push_back(position);
                    }
                }
                if (line[0] == 'v' && line[1] == 'n')
                {
                    std::vector<double> vertexValues = getDoublesOutOfString(line);
                    if (vertexValues.size() == 3)
                    {
                        Vector3D normal = Vector3D(vertexValues[0], vertexValues[1], vertexValues[2]).normalize();
                        normals->push_back(normal);
                    }
                }
                if (line[0] == 'f' && line[1] == ' ')
                {
                    std::vector<int> faceValues = getIntOutOfString(line);
                    if (faceValues.size() == 6)
                    {
                        std::tuple<int, int, int, int, int, int > face = std::tuple<int, int, int, int, int, int>(faceValues[0], faceValues[1], faceValues[2], faceValues[3], faceValues[4], faceValues[5]);
                        faces->push_back(face);
                    }
                    else if (faceValues.size() == 9)
                    {
                        std::tuple<int, int, int, int, int, int > face = std::tuple<int, int, int, int, int, int>(faceValues[0], faceValues[2], faceValues[3], faceValues[5], faceValues[6], faceValues[8]);
                        faces->push_back(face);
                    }
                }
            }
        }
        myfile.close();
        return true;
    }
    return false;
}


void main()
{
    std::vector<Vector3D> positions;
    std::vector<Vector3D> normals;
    std::vector<std::tuple<int, int, int, int, int, int>> faces;

    readOBJ("meshes/cube.obj", &positions, &normals, &faces);

    std::vector<Vector3D> curvedResult;
    int edgeDivisionCount = 100;
    for (std::vector<std::tuple<int, int, int, int, int, int>>::iterator it = faces.begin(); it != faces.end(); ++it)
    {
        std::vector<Vector3D> triangleResult = curvedPNTriangle(edgeDivisionCount, positions[std::get<0>(*it) - 1], positions[std::get<2>(*it) - 1], positions[std::get<4>(*it) - 1], normals[std::get<1>(*it) - 1], normals[std::get<3>(*it) - 1], normals[std::get<5>(*it) - 1]);
        for (std::vector<Vector3D>::iterator itVertex = triangleResult.begin(); itVertex != triangleResult.end(); ++itVertex)
        {
            curvedResult.push_back(*itVertex);
        }
    }
    saveOBJ("meshes/cubeCurved.obj", &curvedResult, nullptr);
}

Okay, so lets see the results:

This discontinuity of the neighboring normal directions of the faces need to be avoided. Now I am not sure, if this is an implementation mistake on my side, or is the curved PN-triangles algorithm the wrong choice to solve my problem. I thought this algorithm would avoid these effects because the results in figure 1 of the paper don't show any of these effects.

In case of correct implementation I would like to ask which algorithm would be a better choice for me?

My requirements are:

1. real time (it must be fast! Calculation per frame)
2. the algorithm must be able to calculate the result within one step (only one shader call need to be enough)
3. supports calculation of normal vectors per vertex (I am not interested in an algorithm which don't deliver the correct normal vector)

Thanks in advance!

Answer 1

As mentioned in the paper PN-triangles are only $C^0$ smooth with respect to adjacent triangles (although they are $G^1$ at vertices). This means that a PN triangle mesh is continuous in position, but are discontinuous in tangent plane and thus in the neighbouring normals on the shared edge of two adjacent triangles. The reason that the surface appears smooth in Figure 1 of the paper is because an additional fake normal field is used (see section 2.1 of the paper). This normal field is continuous and gives the illusion of a smooth surface.

I can point you to a few techniques that might satisfy your requirements

1. Triangular Gregory patches, here the patches will become $G^1$ or tangent plane continous. It able to be implemented in real-time (using for instance tessellation shaders), but it can be quite hard to wrap your head around the Chiyokura-Kimura technique of joining two patches smoothly. The Gregory patch is a triangular B'ezier patch where certain control points are blended between two positions using a rational function. The normal could be analytically determined, but due to the rational functions it becomes hard to take derivatives.
2. Quadratic Approximation of Subdivision surfaces. This technique is able to be implemented in realtime. Although it does not interpolate normal vectors at vertices, it gives smoother results than PN-triangles. It tries to approximate a Loop-subdivision surface with quadratic B'ezier patches.
3. PNG1 triangles Comparable to Gregory patches, but in this case a lot of control points are allowed to shift position. I would advise against implementing this technique as it is very complicated and does not give better results than the standard Gregory patch.
4. Loop subdivision, you can pair PN-triangles with several iterations of Loop subdivision to pre-smooth your mesh. The sharp edges of the triangles will become less and less apparent, but at the cost of increased work to do the subdivision.

Answer 2

Re: do they need additional information about neighboring triangles? If only the edge of the neighbor piece is known you cannot get smooth transitions to the neighbor piece.

Try Polyhedral Splines