A short example on how triangle normals can be used to compute vertex normals for either smooth or flat shading.

triangle_to_vertex_normals.cpp

// Please note: The code below has not been compiled. Includes are missing.
// This snippet is just for illustration purposes. Feedback is highly appreciated, of course.

// First we'll read the stl data into local arrays. Each triplet in `rawPoints` corresponds to one
// 3D point and each triplet in `rawTriNormals` to one 3D triangle normal.
std::vector <float> rawPoints, rawTriNormals;
std::vector <unsigned int> triIndices, solids;

stl_reader::ReadStlFile ("geometry.stl", rawPoints, rawTriNormals, triIndices, solids);

const size_t numTris = triIndices.size() / 3;

// We'll use vector types from the glm library for more readable math.
std::vector <glm::vec3> points;
for (size_t i = 0; i + 2 < rawPoints.size (); i += 3)
  points.push_back ({rawPoints [i], rawPoints [i+1], rawPoints [i+2]});

std::vector <glm::vec3> triNormals;
for (size_t i = 0; i + 2 < rawTriNormals.size (); i += 3)
  triNormals.push_back ({rawTriNormals [i], rawTriNormals [i+1], rawTriNormals [i+2]});

// In `smoothVertexNormals` we'll store the normal of each vertex. They are computed
// by first summing up the normals of adjacent triangles and then normalizing the
// resulting normal to get a unit normal for each vertex.
// Those `smoothVertexNormals` can then be used alongside `points` and `triIndices` to draw indexed
// primitives using `glDrawElements` or similar methods.
// While smooth vertex normals lead to a smooth and round appearance of objects, they are are not
// well suited to represent the surface orientation at crease vertices, i.e., vertices at which
// triangles meet at large angles. See `flatPoints` and `flatNormals` below for an alternate approach.
std::vector <glm::vec3> smoothVertexNormals (points.size (), glm::vec3 (0));
for (size_t itri = 0; itri < numTris; ++itri)
{
  const size_t baseIndex = itri * 3;
  for (size_t icorner = 0; icorner < 3; ++icorner)
    smoothVertexNormals [triIndices [baseIndex + icorner]] += triNormals [itri];
}

for (auto& n : triNormals)
  n = glm::normalize (n);

// Sometimes flat shading or partial flat shading is preferrable to smooth shading. While OpenGL
// offers a shading mode for that, the results from that mode are often not accurate, since only
// the normal of the first vertex of each triangle is considered for lighting calculations. Here is
// a simple example on how to set up points and normals for accurate flat shading. Note that accurate
// flat shading can also be achieved using geometry shaders. But this is probably out of the scope
// of this example...
// To render the resulting `flatPoints` and `flatNormals`, the `triIndices` array can now be ignored
// and a method like `glDrawArrays` should be used.
std::vector <glm::vec3> flatPoints, flatNormals;
for (size_t itri = 0; itri < numTris; ++itri)
{
  const size_t baseIndex = itri * 3;
  for (size_t icorner = 0; icorner < 3; ++icorner)
  {
    flatPoints.push_back (points [triIndices [baseIndex + icorner]]);
    flatNormals.push_back (triNormals [itri]);
  }
}

This works well but only after commenting out RemoveDoubles and adding coordsOut.push_back(c[i]); in the loop that reads the position.

With RemoveDoubles in code

After commenting out RemoveDoubles and adding coordsOut.push_back(c[i]); in the loop that reads the position