How to unify normal orientation

Question

I've been trying to realize a mesh that has all face normals pointing outward. In order to realize this, I load a mesh from a *.ctm file, then walk over all triangles to determine the normal using a cross product and if the normal is pointing to the negative z direction, I flip v1 and v2 (thus the normal orientation). After this is done I save the result to a *.ctm file and view it with Meshlab.

The result in Meshlab still shows that normals are pointing in both positive and negative z direction ( can be seen from the black triangles). Also when viewing the normals in Meshlab they are really pointing backwards.

Can anyone give me some advice on how to solve this?

The source code for the normalization part is:

pcl::PointCloud<pcl::PointXYZRGBA>::Ptr cloud1 (new pcl::PointCloud<pcl::PointXYZRGBA> ());
pcl::fromROSMsg (meshFixed.cloud,*cloud1);for(std::vector<pcl::Vertices>::iterator it = meshFixed.polygons.begin(); it != meshFixed.polygons.end(); ++it)
{
    alglib::real_2d_array v0;
    double _v0[] = {cloud1->points[it->vertices[0]].x,cloud1->points[it->vertices[0]].y,cloud1->points[it->vertices[0]].z};
    v0.setcontent(3,1,_v0); //3 rows, 1col
    alglib::real_2d_array v1;
    double _v1[] = {cloud1->points[it->vertices[1]].x,cloud1->points[it->vertices[1]].y,cloud1->points[it->vertices[1]].z};
    v1.setcontent(3,1,_v1); //3 rows, 1col
    alglib::real_2d_array v2;
    double _v2[] = {cloud1->points[it->vertices[2]].x,cloud1->points[it->vertices[2]].y,cloud1->points[it->vertices[2]].z};
    v2.setcontent(1,3,_v2); //3 rows, 1col
    alglib::real_2d_array normal;
    normal = cross(v1-v0,v2-v0);
    //if z<0 change indices order v1->v2 and v2->v1
    alglib::real_2d_array normalizedNormal;
    if(normal[2][0]<0)
    {
            int index1,index2;
            index1 = it->vertices[1];
            index2 = it->vertices[2];
            it->vertices[1] = index2;
            it->vertices[2] = index1;
            //make normal of length 1
            double normalScaling = 1.0/sqrt(dot(normal,normal));
            normal[0][0] = -1*normal[0][0];
            normal[1][0] = -1*normal[1][0];
            normal[2][0] = -1*normal[2][0];
            normalizedNormal = normalScaling * normal;
    }
    else
    {
            //make normal of length 1
            double normalScaling = 1.0/sqrt(dot(normal,normal));
            normalizedNormal = normalScaling * normal;
    }
    //add to normal cloud
    pcl::Normal pclNormalizedNormal;
    pclNormalizedNormal.normal_x = normalizedNormal[0][0];
    pclNormalizedNormal.normal_y = normalizedNormal[1][0];
    pclNormalizedNormal.normal_z = normalizedNormal[2][0];
    normalsFixed.push_back(pclNormalizedNormal);
} 

The result from this code is:

I've found some code in the VCG library to orient the face and vertex normals. After using this a large part of the mesh has correct face normals, but not all.

The new code:

// VCG library implementation
    MyMesh m;
    // Convert pcl::PolygonMesh to VCG MyMesh
    m.Clear();
    // Create temporary cloud in to have handy struct object
    pcl::PointCloud<pcl::PointXYZRGBA>::Ptr cloud1 (new pcl::PointCloud<pcl::PointXYZRGBA> ());
    pcl::fromROSMsg (meshFixed.cloud,*cloud1);
    // Now convert the vertices to VCG MyMesh
    int vertCount = cloud1->width*cloud1->height;
    vcg::tri::Allocator<MyMesh>::AddVertices(m, vertCount);
    for(unsigned int i=0;i<vertCount;++i)
        m.vert[i].P()=vcg::Point3f(cloud1->points[i].x,cloud1->points[i].y,cloud1->points[i].z);
    // Now convert the polygon indices to VCG MyMesh => make VCG faces..
    int triCount = meshFixed.polygons.size();
    if(triCount==1)
    {
        if(meshFixed.polygons[0].vertices[0]==0 && meshFixed.polygons[0].vertices[1]==0 && meshFixed.polygons[0].vertices[2]==0)
            triCount=0;
    }
    Allocator<MyMesh>::AddFaces(m, triCount);
    for(unsigned int i=0;i<triCount;++i)
    {
        m.face[i].V(0)=&m.vert[meshFixed.polygons[i].vertices[0]];
        m.face[i].V(1)=&m.vert[meshFixed.polygons[i].vertices[1]];
        m.face[i].V(2)=&m.vert[meshFixed.polygons[i].vertices[2]];
    }

    vcg::tri::UpdateBounding<MyMesh>::Box(m);
    vcg::tri::UpdateNormal<MyMesh>::PerFace(m);
    vcg::tri::UpdateNormal<MyMesh>::PerVertexNormalizedPerFace(m);
    printf("Input mesh  vn:%i fn:%i\n",m.VN(),m.FN());

    // Start to flip all normals to outside
    vcg::face::FFAdj<MyMesh>::FFAdj();
    vcg::tri::UpdateTopology<MyMesh>::FaceFace(m);
    bool oriented, orientable;
    if ( vcg::tri::Clean<MyMesh>::CountNonManifoldEdgeFF(m)>0 ) {
        std::cout << "Mesh has some not 2-manifold faces, Orientability requires manifoldness" << std::endl; // text
        return; // can't continue, mesh can't be processed
    }
    vcg::tri::Clean<MyMesh>::OrientCoherentlyMesh(m, oriented,orientable);
    vcg::tri::Clean<MyMesh>::FlipNormalOutside(m);
    vcg::tri::Clean<MyMesh>::FlipMesh(m);
    //vcg::tri::UpdateTopology<MyMesh>::FaceFace(m);
    //vcg::tri::UpdateTopology<MyMesh>::TestFaceFace(m);
    vcg::tri::UpdateNormal<MyMesh>::PerVertexNormalizedPerFace(m);
    vcg::tri::UpdateNormal<MyMesh>::PerVertexFromCurrentFaceNormal(m);

    // now convert VCG back to pcl::PolygonMesh
    pcl::PointCloud<pcl::PointXYZRGBA>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZRGBA>);
    cloud->is_dense = false;
    cloud->width = vertCount;
    cloud->height = 1;
    cloud->points.resize (vertCount);
    // Now fill the pointcloud of the mesh
    for(int i=0; i<vertCount; i++)
    {
        cloud->points[i].x = m.vert[i].P()[0];
        cloud->points[i].y = m.vert[i].P()[1];
        cloud->points[i].z = m.vert[i].P()[2];
    }
    pcl::toROSMsg(*cloud,meshFixed.cloud);
    std::vector<pcl::Vertices> polygons;
    // Now fill the indices of the triangles/faces of the mesh
    for(int i=0; i<triCount; i++)
    {
        pcl::Vertices vertices;
        vertices.vertices.push_back(m.face[i].V(0)-&*m.vert.begin());
        vertices.vertices.push_back(m.face[i].V(1)-&*m.vert.begin());
        vertices.vertices.push_back(m.face[i].V(2)-&*m.vert.begin());
        polygons.push_back(vertices);
    }
    meshFixed.polygons = polygons;

Which results in: (Meshlab still shows normals are facing both sides)

Answer 1

I finally solved the problem. So I'm still using VCG library. From the above new code I slightly updated the following section:

vcg::tri::Clean<MyMesh>::OrientCoherentlyMesh(m, oriented,orientable);
//vcg::tri::Clean<MyMesh>::FlipNormalOutside(m);
//vcg::tri::Clean<MyMesh>::FlipMesh(m);
//vcg::tri::UpdateTopology<MyMesh>::FaceFace(m);
//vcg::tri::UpdateTopology<MyMesh>::TestFaceFace(m);
vcg::tri::UpdateNormal<MyMesh>::PerVertexNormalizedPerFace(m);
vcg::tri::UpdateNormal<MyMesh>::PerVertexFromCurrentFaceNormal(m);

Now I've updated the vcg::tri::Clean<MyMesh>::OrientCoherentlyMesh() function in clean.h. Here the update is to orient the first polygon of a group correctly. Also after swapping the edge the normal of the face is calculated and updated.

static void OrientCoherentlyMesh(MeshType &m, bool &Oriented, bool &Orientable)
{
    RequireFFAdjacency(m);
    assert(&Oriented != &Orientable);
    assert(m.face.back().FFp(0));    // This algorithms require FF topology initialized

    Orientable = true;
    Oriented = true;

    tri::UpdateSelection<MeshType>::FaceClear(m);
    std::stack<FacePointer> faces;

    for (FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
    {
        if (!fi->IsD() && !fi->IsS())
        {
            // each face put in the stack is selected (and oriented)
            fi->SetS();
            // New section of code to orient the initial face correctly
            if(fi->N()[2]>0.0)
            {
                face::SwapEdge<FaceType,true>(*fi, 0);
                face::ComputeNormal(*fi);
            }
            // End of new code section.
            faces.push(&(*fi));

            // empty the stack
            while (!faces.empty())
            {
                FacePointer fp = faces.top();
                faces.pop();

                // make consistently oriented the adjacent faces
                for (int j = 0; j < 3; j++)
                {
                   //get one of the adjacent face
                   FacePointer fpaux = fp->FFp(j);
                   int iaux = fp->FFi(j);

                   if (!fpaux->IsD() && fpaux != fp && face::IsManifold<FaceType>(*fp, j))
                   {              
                      if (!CheckOrientation(*fpaux, iaux))
                      {
                          Oriented = false;

                          if (!fpaux->IsS())
                          {
                               face::SwapEdge<FaceType,true>(*fpaux, iaux);
                               // New line to update face normal
                               face::ComputeNormal(*fpaux);
                               // end of new section.
                               assert(CheckOrientation(*fpaux, iaux));
                          }
                          else
                          {
                               Orientable = false;
                               break;
                          }
                       }

                       // put the oriented face into the stack

                       if (!fpaux->IsS())
                       {
                            fpaux->SetS();
                            faces.push(fpaux);
                       }
                   }
               }
           }
       }
       if (!Orientable) break;
    }
}

Besides I also updated the function bool CheckOrientation(FaceType &f, int z) to perform a calculation based on normal z-direction.

template <class FaceType>
bool CheckOrientation(FaceType &f, int z)
{
    // Added next section to calculate the difference between normal z-directions
    FaceType *original = f.FFp(z);
    double nf2,ng2;
    nf2=f.N()[2];
    ng2=original->N()[2];
    // End of additional section
    if (IsBorder(f, z))
        return true;
    else
    {
        FaceType *g = f.FFp(z);
        int gi = f.FFi(z);
        // changed if statement from: if (f.V0(z) == g->V1(gi))
        if (nf2/abs(nf2)==ng2/abs(ng2))
            return true;
        else
            return false;
    }
}

The result is as I expect and desire from the algorithm: