random points in a hexarea / shape

Question

I have a simple hexagonal grid, where I'm selecting a group of hexagons, which will then be filled with some random points.

Let me explain the exact procedure of generating the points:

• Using a list of hex coordinates I'm selecting hexagons.
• Group hexagons into areas.
• Store all edge points, for each area individually.
• Loop through areas:
  • Calculate area boundary ( required for definig the range for random point generator )
  • Draw area (bounding) square. (to see whether the computation was correct)
  • Based on the bounding square min,max values generate random points
  • Test if point is inside the area shape. (The shape is defined by area edge points)
  • If the point passed the above test, it's the pushed into an array.
• loop through the array of points and draw them on screen.

Now I tried these two methods to determine whether a point lies inside a specific shape.

cn_PnPoly: function( P, V, n ){ ///// Point , array of vertices , array size ////
        var cn = 0, vt = 0.0;
        for (var i=0; i< (n-1); i++) {    // edge from V[i]  to V[i+1]
           if (((V[i].y <= P.y) && (V[i+1].y > P.y))     // an upward crossing
            || ((V[i].y > P.y) && (V[i+1].y <=  P.y))) { // a downward crossing
                // compute  the actual edge-ray intersect x-coordinate
                vt = (P.y  - V[i].y) / (V[i+1].y - V[i].y);
                if (P.x <  V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect
                     ++cn;   // a valid crossing of y=P.y right of P.x
            }
        }
        return (cn&1);    // 0 if even (out), and 1 if  odd (in)
    },

result:

isInHexBounary: function( p, points, l ){ ///// Point , array of vertices , array size ////
        var result = false;
          for (i = 0, j = l - 1; i < l; j = i++) {
            if ((points[i].y > p.y) != (points[j].y > p.y) && (p.x < (points[j].x - points[i].x) * (p.y - points[i].y) / (points[j].y-points[i].y) + points[i].x)) {
                result = !result;
             }
          }
          return result;
    },

result:

I suppose the first method requires to have all the points in specific order, and that's why it's not working properly. But the second one seems to be working almost correct , apart from some parts. Any idea what I'm doing wrong?

Update:

It's turning out that for mos tof the algorhytms it's required to have the points in a specific order, so I calculated angle of each point relative to the average center point. And then sorted all the points by angle. Both algorithms now return similar results. Although there is sitll some points leaking through the area shape.

findCenter: function( points ){
        var x = 0, y = 0, i, len = points.length;
        for (i = 0; i < len; i++) {
            x += points[i].x;
            y += points[i].y;
        }
        return {x: x / len, y: y / len};
    },

    findAngles: function(c, points){
        var i, len = points.length, p, dx, dy;
        for (i = 0; i < len; i++) {
            p = points[i];
            dx = p.x - c.x;
            dy = p.y - c.y;
            p.angle = Math.atan2(dy, dx);
        }
    },

    sortByAngle: function( points ){
        points.sort(function(a, b) {
          if (a.angle > b.angle) return 1;
          else if (a.angle < b.angle) return -1;
          return 0;
        });
    },

result:

Answer 1

Frankly, I would make it simpler.

1. Make sampling in one single hexagon. Just make bounding box with size (2s,2s cos(30)), sample (x,y) and reject it if it is outside the hexagon. Efficiency should be 3/4 (check hexagon area). Position this sampling hexagon at (0,0) such that sampling is very easy, and, what's more important, very TESTABLE

2. For each area maintain array of hexagon centers, say, with size N

3. Sample in two steps. First, sample integer 1...N and select what hexagon in the area is about to be selected. Second, sample from your (x,y) from step #1, from your single hexagon at (0,0). Last, shift sampled (x,y) by randomly selected hexagon center

UPDATE

Actually, I believe there is a way to sample point (x,y) in a single hexagon with 100% efficiency without any rejection/acceptance. Look at picture above. Red is whole hexagon, and rectangular blue area is where you sample points. If sampled point is within red area you take it and move on. If it is not in the red and inside blue A triangle, you map it into black A' triangle. If point is not in the red but in the blue B triangle, you remap it into black B' triangle

Remapping is quite simple linear transformation.

Finally you have sampling with three random numbers as input (one used to select destination hexagon, and two used to sample random point) and it will guaranteed to return random point somewhere in the desired area

UPDATE II

As noted by Denis Sheremet, better mapping would be A->B' and B->A'. Assuming hexagon center at (0,0), overall just two reflections - one over center and another over middle of the triangle

Answer 2

I implemented Severin Pappadeux's awesome answer in C++ (but should be trivially convertible to any other language):

    // randA and randB must be in the range [0, 1]
    static glm::vec2 SampleHexagon(float randA, float randB) {
        // Algorithm based on https://stackoverflow.com/a/39262805/3841944 by user Severin Pappadeux
        //
        // Hexagon map:
        //  ___________________
        // |   /|         |\   |
        // |  / |         | \b'|
        // | / a|         |A \ |
        // |/___|    C    |___\|
        // |\   |         |   /|
        // | \ b|         |B / |
        // |  \ |         | /a'|
        // |   \___________/___|
        // |    |         | <- "startA"
        // |    |         |<-->| <- "widthA"
        // |    | <- "startC"
        // |<----------------->| "width"
        //
        // C = center part
        // a, b = left (upper and lower) part
        // A, B = right (upper and lower) part
        // a' -> area that will be remapped to a
        // b' -> area that will be remapped to b

        // The algorithm is to:
        // 1. Generate random points in the rectangle spanning C, A, B, b' and a'
        // 2. Move all points in a' and b' into a and b
        // 3. (Optional) Make hexagon regular
        // 4. (Optional) Remap coordinates to range [-1, 1]

        // Coordinates are in the range [0, 1]
        const float width  = 1;
        const float widthC = width / 2;
        const float widthA = (width - widthC) / 2;
        const float startC = widthA;
        const float startA = startC + widthC;
        const float slope  = .5f / widthA;


        // 1. Remap x into the rectangle spanning C, A, B, b' and a'
        float x = startC + randA * .75f;
        float y = randB;


        // 2. Move all points in a' and b' into a and b

        // If we are in the last third (A, b', B and a')
        if (startA < x) {
            float localX = x - startA;

            if (y > .5f) {  // And the upper half of it (in A or b')
                float localY = y - .5f;
                if (localY > .5f - localX * slope) { // And above the diagonal
                    // Move the point so it is in the triangle to 'b'
                    x -= startA;
                    y -= .5f;
                }
            } else { // And the lower half of it (in B or a')
                float localY = y;
                if (localY < localX * slope) {  // And we are below the diagonal
                    // Move the point so it is in the triangle to 'a'
                    x -= startA;
                    y += .5f;
                }
            }
        }

        // 3. Make to regular hexagon (currently the hexagon is too high, because we assumed width == height)

        // Divide the hexagon into 6 regular triangles and use the Pythagorean theorem, and this calculation should be fairly obvious
        float regularHexagonWidthToHeightRatio = std::sqrt(1.f*1.f - .5f*.5f);
        y = (y - .5f) * regularHexagonWidthToHeightRatio + .5f; // Center around 0, scale, center to .5

        // 4. Remap coordinates to -1 to 1 (optional)
        x = x * 2 - 1;
        y = y * 2 - 1;

        return { x, y };
    }

Again many thanks to Severin Pappadeux, upvote his answer.