Procedural generation of stars with skybox

I am attempting to procedurally generate a star-filled background in OpenGL.

The approach I am taking is to create a skybox with a cubemap texture. Each side of the cubemap texture essentially consists of a 2048x2048 black image with randomly selected texels set to White. Here is the result:

I'm not sure how obvious it is from the image, but when moving around a very distinct box shape can be made out as stars close to the edge of the box appear smaller and closer together. How can I prevent this? Do I need to abandon the skybox approach and use something like a skysphere instead?

EDIT: here is how I am mapping the cubemap onto the sky.

// Create and bind texture.
glGenTextures(1, &texture_);
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_CUBE_MAP, texture_);

for (unsigned int i = 0; i < 6; ++i) {
    std::vector<std::uint8_t> image = generateTexture(TEXTURE_WIDTH, TEXTURE_HEIGHT);
    glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, 0, GL_RGB, TEXTURE_WIDTH, TEXTURE_HEIGHT,
                 0, GL_RGB, GL_UNSIGNED_BYTE, image.data());
}

// Set texture parameters.
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_CUBE_MAP, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);

Here is the definition of the generateTexture function:

std::vector<std::uint8_t> Stars::generateTexture(GLsizei width, GLsizei height) {
    std::vector<std::uint8_t> image(static_cast<std::size_t>(3 * width * height));

    add_stars(image, NUM_STARS);

    return image;
}

void Stars::add_stars(std::vector<std::uint8_t>& image, unsigned int nStars) {
    std::default_random_engine eng;
    std::uniform_int_distribution<std::size_t> dist(0, image.size() / 3 - 1);

    while (nStars--) {
        std::size_t index = 3 * dist(eng);

        image[index++] = 255;
        image[index++] = 255;
        image[index++] = 255;
    }
}

EDIT2: here is the draw function used to render the sky.

void Stars::draw(const Camera& camera) const {
    // Skybox will be rendered last. In order to ensure that the stars are rendered at the back of
    // the scene, the depth buffer is filled with values of 1.0 for the skybox -- this is done in
    // the vertex shader. We need to make sure that the skybox passes the depth te3t with values
    // less that or equal to the depth buffer.
    glDepthFunc(GL_LEQUAL);

    program_.enable();

    // Calculate view-projection matrix and set the corresponding uniform. The view matrix must be
    // stripped of translation components so that the skybox follows the camera.
    glm::mat4 view       = glm::mat4(glm::mat3(camera.viewMatrix()));
    glm::mat4 projection = camera.projectionMatrix();

    glm::mat4 VP = projection * view;
    glUniformMatrix4fv(program_.uniformLocation("VP"), 1, GL_FALSE, glm::value_ptr(VP));

    // Bind buffer objects and texture to current context and draw.
    glBindVertexArray(vao_);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo_);
    glBindTexture(GL_TEXTURE_CUBE_MAP, texture_);

    glDrawElements(GL_TRIANGLES, static_cast<GLsizei>(INDICES.size()), GL_UNSIGNED_INT,
                   reinterpret_cast<GLvoid *>(0));

    glBindVertexArray(0);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
    glBindTexture(GL_TEXTURE_CUBE_MAP, 0);
    program_.disable();

    glDepthFunc(GL_LESS);
}

generate stars uniformly in some cubic volume

x=2.0*Random()-1.0; // <-1,+1>
y=2.0*Random()-1.0; // <-1,+1>
z=2.0*Random()-1.0; // <-1,+1>

project them on unit sphere

So just compute the length of vector (x,y,z) and divide the coordinates by it.
project the result onto the cube map

Each side of cube is defined by the plane so find intersection of ray casted from (0,0,0) through Cartesian star position and the planes. Take the intersection with shortest distance to (0,0,0) and use that as final star position.

The implementation could be something like this OpenGL&C++ code:

    glClearColor(0.0,0.0,0.0,0.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    int i,n=10000;
    float a,b,x,y,z;
    //RandSeed=8123456789;
    n=obj.pnt.num;  // triangulated sphere point list

    glDepthFunc(GL_LEQUAL);
    glEnable(GL_BLEND);
    glBlendFunc(GL_ONE,GL_ONE);

    glPointSize(2.0);
    glBegin(GL_POINTS);
    for (i=0;i<n;i++)
        {
        // equidistant points instead of random to test this
        x=obj.pnt[i].p[0];
        y=obj.pnt[i].p[1];
        z=obj.pnt[i].p[2];
/*
        // random star spherical position
        a=2.0*M_PI*Random();
        b=M_PI*(Random()-0.5);
        // spherical 2 cartessian r=1;
        x=cos(a)*cos(b);
        y=sin(a)*cos(b);
        z=       sin(b);
*/
        // redish sphere map
        glColor3f(0.6,0.3,0.0); glVertex3f(x,y,z);
        // cube half size=1 undistort // using similarities like: yy/xx = y/x
             if ((fabs(x)>=fabs(y))&&(fabs(x)>=fabs(z))){ y/=x; z/=x; if (x>=0) x=1.0; else x=-1.0; }
        else if ((fabs(y)>=fabs(x))&&(fabs(y)>=fabs(z))){ x/=y; z/=y; if (y>=0) y=1.0; else y=-1.0; }
        else if ((fabs(z)>=fabs(x))&&(fabs(z)>=fabs(y))){ x/=z; y/=z; if (z>=0) z=1.0; else z=-1.0; }
        // bluish cube map
        glColor3f(0.0,0.3,0.6); glVertex3f(x,y,z);
        }
    glEnd();
    glPointSize(1.0);
    glDisable(GL_BLEND);
    glFlush();
    SwapBuffers(hdc);

Looks like it works as it should here preview (of the blended sphere/cube map):

Although it looks like there are holes but there are none (it is may be some blend error) if I disable the sphere map render then there are no visible holes or distortions in the mapping.

The sphere triangulation mesh obj used to test this is taken from here:

Sphere triangulation

[Edit1] yes there was a silly blending error

I repaired the code ... but the problem persists anyway. does not matter this mapping is working as should here the updated code result:

So just adapt the code to your texture generator ...

[Edit2] Random stars

glClearColor(0.0,0.0,0.0,0.0);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
int i;
float x,y,z,d;
RandSeed=8123456789;

glDepthFunc(GL_LEQUAL);
glEnable(GL_BLEND);
glBlendFunc(GL_ONE,GL_ONE);

glPointSize(2.0);
glBegin(GL_POINTS);
for (i=0;i<1000;i++)
    {
    // uniform random cartesian stars inside cube
    x=(2.0*Random())-1.0;
    y=(2.0*Random())-1.0;
    z=(2.0*Random())-1.0;
    // project on unit sphere
    d=sqrt((x*x)+(y*y)+(z*z));
    if (d<1e-3) { i--; continue; }
    d=1.0/d;
    x*=d; y*=d; z*=d;
    // redish sphere map
    glColor3f(0.6,0.3,0.0); glVertex3f(x,y,z);
    // cube half size=1 undistort using similarities like: y/x = y'/x'
         if ((fabs(x)>=fabs(y))&&(fabs(x)>=fabs(z))){ y/=x; z/=x; if (x>=0) x=1.0; else x=-1.0; }
    else if ((fabs(y)>=fabs(x))&&(fabs(y)>=fabs(z))){ x/=y; z/=y; if (y>=0) y=1.0; else y=-1.0; }
    else if ((fabs(z)>=fabs(x))&&(fabs(z)>=fabs(y))){ x/=z; y/=z; if (z>=0) z=1.0; else z=-1.0; }
    // bluish cube map
    glColor3f(0.0,0.3,0.6); glVertex3f(x,y,z);
    }
glEnd();
glPointSize(1.0);
glDisable(GL_BLEND);
glFlush();
SwapBuffers(hdc);

Here Blend of booth (1000 stars):

And Here only the cube-map (10000 stars)

[Edit3] The Blend problem solved

It was caused by Z-fighting and occasional changing of sign for some coordinates during the projection due to forgotten fabs here fixed code:

    glClearColor(0.0,0.0,0.0,0.0);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    int i;
    float x,y,z,d;
    RandSeed=8123456789;

    glDepthFunc(GL_ALWAYS);
//  glDepthFunc(GL_LEQUAL);
    glEnable(GL_BLEND);
    glBlendFunc(GL_ONE,GL_ONE);

    glPointSize(2.0);
    glBegin(GL_POINTS);
    for (i=0;i<25000;i++)
        {
        // uniform random cartesian stars inside cube
        x=(2.0*Random())-1.0;
        y=(2.0*Random())-1.0;
        z=(2.0*Random())-1.0;
        // project on unit sphere
        d=sqrt((x*x)+(y*y)+(z*z));
        if (d<1e-3) { i--; continue; }
        d=1.0/d;
        x*=d; y*=d; z*=d;
        // redish sphere map
        glColor3f(0.6,0.3,0.0); glVertex3f(x,y,z);
        // cube half size=1 undistort using similarities like: y/x = y'/x'
             if ((fabs(x)>=fabs(y))&&(fabs(x)>=fabs(z))){ y/=fabs(x); z/=fabs(x); if (x>=0) x=1.0; else x=-1.0; }
        else if ((fabs(y)>=fabs(x))&&(fabs(y)>=fabs(z))){ x/=fabs(y); z/=fabs(y); if (y>=0) y=1.0; else y=-1.0; }
        else if ((fabs(z)>=fabs(x))&&(fabs(z)>=fabs(y))){ x/=fabs(z); y/=fabs(z); if (z>=0) z=1.0; else z=-1.0; }
        // bluish cube map
        glColor3f(0.0,0.3,0.6); glVertex3f(x,y,z);
        }
    glEnd();
    glPointSize(1.0);
    glDisable(GL_BLEND);
    glFlush();
    SwapBuffers(hdc);

And here the Blend result finally the colors are as should be so the sphere and cube stars overlaps perfectly (white) while viewing from (0,0,0):