Including edges in this pnpoly algorithm

Question

I have this point in polygon function to use in my pathfinding program.

int point_in_pol(int vertcount, float *vertx, float *verty, int vertexx, int vertexy){
    double vertexx1;
    vertexx1 = vertexx;
    double vertexy1;
    vertexy1 = vertexy;
  int i ,j, c = 0;
  for (i = 0, j = vertcount-1; i < vertcount; j = i++) {
    if ( (((verty[i]>=vertexy1) && (verty[j]<=vertexy1) )  ||  ((verty[i]<=vertexy1)   && (verty[j]>=vertexy1) )) &&
     (vertexx1 < (vertx[j]-vertx[i]) * (vertexy1-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

This function returns true if the point is in the polygon. However, it does not behave properly if the point given is on the edge of the polygon. What changes should i do here to make it return true if the point is in the edge?

whole code:

#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>


typedef struct Node Node;
typedef struct Qnode Qnode;
void enqueue(Node* node);
void enqueue_left(Node* node);
Node* generate(int x, int y);
Node* dequeue();
void expand(Node* node, int xend, int yend);
int point_in_pol(int vertcount, float *vertx, float *verty, int vertexx, int vertexy);
struct Node{
    Node* predecessor;
    Node* up;
    Node* down;
    Node* left;
    Node* right;
    int xpos;
    int ypos;
    int marked;
};
struct Qnode{
    Qnode* next;
    Node* Gnode; 
};
Qnode* front = NULL;
Qnode* rear = NULL;
int queuesize = 0;
int des;
int solutioncost = 0;
Node* closednodes[80000];
int nodesclosed = 0;
float polygonx[20][50];
float polygony[20][50];
int polycount = 0, vertcount = 0;
int vertcounts[20];

void enqueue(Node* node){
    if (queuesize == 0){
        rear = (Qnode*)malloc(sizeof(Qnode));
        rear->Gnode = node;
        rear->next = NULL;
        front = rear;
    }
    else{
        Qnode* temp = (Qnode*)malloc(sizeof(Qnode));
        rear->next = temp;
        temp->Gnode = node;
        temp->next = NULL;
        rear = temp;
    }
        queuesize++;
}
void enqueue_left(Node* node){
    if(queuesize == 0){
        front = (Qnode*)malloc(sizeof(Qnode));
        front->Gnode = node;
        front->next = NULL;
        rear = front;
    }
    else{
        Qnode* temp = (Qnode*)malloc(sizeof(Qnode));
        temp->Gnode = node;
        temp->next = front;
        front = temp;
    }
    queuesize++;
}

Node* dequeue(){
    Qnode* temp = front;
    if (queuesize == 0){
        printf("Error!");
    }
    else{
        Node* temp1 = front->Gnode;
        temp = front->next;
        free(front);
        front = temp;
        queuesize--;
        return temp1;
    }

}

Node* generate(int x, int y){
    int i = 0, j = 0;
    //printf("Generating node (%d,%d)...\n", x, y);
    if ((x >0 && x <=400) && (y >0 && y <=200)){
    Node* temp = (Node*)malloc(sizeof(Node));
    temp->xpos = x;
    temp->ypos = y;
    temp->marked = 0;
    for(i; i<polycount; i++){
        if(point_in_pol(vertcounts[i], polygonx[i],polygony[i], x, y)){
            temp->marked = 1;
        }
    }
    temp->up = NULL;
    temp->predecessor = NULL;
    temp->down = NULL;
    temp->left = NULL;
    temp->right = NULL;
    return temp;
    }
    else{
        printf("Invalid starting point! \n");
    }
}

void expand(Node* node, int xend, int yend){
    //printf("Expanding node (%d, %d)...\n", node->xpos, node->ypos);
    solutioncost++;
    int flag = 0;
    closednodes[nodesclosed] = node;
    nodesclosed++;
    dequeue();
    if(node->marked == 1){
    //  printf("Cannot expand. Part of a polygon.\n");
    }

    else{
        if (node->xpos == xend && node->ypos == yend){
            printf("Node reached!");
            printf("Path in reverse order: \n(%d, %d)\n", xend, yend);
            while(node->predecessor!= NULL){
                printf("(%d, %d)\n", node->predecessor->xpos, node->predecessor->ypos);
                node = node->predecessor;
            }
            queuesize = 0; 
            return;
        }
        if (node->xpos-1 >0 && node->left == NULL){
            int k = 0;
            int j = 0;
            Qnode* temp2 = front;
            for(k; k<queuesize; k++){
                int xx = temp2->Gnode->xpos;
                int yy = temp2->Gnode->ypos;
                if (xx == node->xpos-1 && yy == node->ypos)
                    flag = 1;
                temp2 = temp2->next;
                }
            for(j; j<nodesclosed; j++){
                int xx = closednodes[j]->xpos;
                int yy = closednodes[j]->ypos;
                if (xx == node->xpos-1 && yy == node->ypos)
                    flag = 1;
            }
            if (flag == 0){
                node->left = generate(node->xpos-1, node->ypos);
                node->left->predecessor = node;
                node->left->right = node;
                if(des == 1)
                    enqueue(node->left);
                else if(des == 2)
                    enqueue_left(node->left);
            }
            else{
                //printf("(%d, %d) not generated.\n", node->xpos-1, node->ypos);
                flag = 0;
            }
        }
        if (node->xpos+1 <=400 && node->right ==NULL){
        int k = 0;
        int j = 0;
        Qnode* temp2 = front;
            for(k; k<queuesize; k++){
                int xx = temp2->Gnode->xpos;
                int yy = temp2->Gnode->ypos;
                if (xx == node->xpos+1 && yy == node->ypos)
                    flag = 1;
                temp2 = temp2->next;
                }

            for(j; j<nodesclosed; j++){
                int xx = closednodes[j]->xpos;
                int yy = closednodes[j]->ypos;
                if (xx == node->xpos+1 && yy == node->ypos)
                    flag = 1;
            }
            if (flag == 0){
                node->right = generate(node->xpos+1, node->ypos);
                node->right->predecessor = node;
                node->right->left = node;
                if(des == 1)
                    enqueue(node->right);
                else if(des == 2)
                    enqueue_left(node->right);
            }
            else{
                //printf("(%d, %d) not generated.\n", node->xpos+1, node->ypos);
                flag = 0;
            }
        }
        if (node->ypos+1 <=200 && node->up ==NULL){
        int k = 0;
        int j = 0;
        Qnode* temp2 = front;
        for(k; k<queuesize; k++){
            int xx = temp2->Gnode->xpos;
            int yy = temp2->Gnode->ypos;
            if (xx == node->xpos && yy == node->ypos+1)
                flag = 1;
            temp2 = temp2->next;
                }   
            for(j; j<nodesclosed; j++){
                int xx = closednodes[j]->xpos;
                int yy = closednodes[j]->ypos;
                if (xx == node->xpos && yy == node->ypos+1)
                    flag = 1;
            }

            if (flag ==0){
                node->up = generate(node->xpos, node->ypos+1);
                node->up->predecessor = node;
                node->up->down = node;
            if(des == 1)
                enqueue(node->up);
            else if(des == 2)
                enqueue_left(node->up);
            }
            else{
                //printf("(%d, %d) not generated.\n", node->xpos, node->ypos+1);
            flag = 0;
            }
        }
        if (node->ypos-1 >0 && node->down ==NULL){
        int k = 0;
        int j = 0;
        Qnode* temp2 = front;
            for(k; k<queuesize; k++){
                int xx = temp2->Gnode->xpos;
                int yy = temp2->Gnode->ypos;
                if (xx == node->xpos && yy == node->ypos-1)
                    flag = 1;
                temp2 = temp2->next;
                }

            for(j; j<nodesclosed; j++){
                int xx = closednodes[j]->xpos;
                int yy = closednodes[j]->ypos;
                if (xx == node->xpos && yy == node->ypos-1)
                    flag = 1;
            }
            if (flag ==0){
                node->down = generate(node->xpos, node->ypos-1);
                node->down->predecessor = node;
                node->down->up = node;
            if(des == 1)
                enqueue(node->down);
            else if(des == 2)
                enqueue_left(node->down);
            }
            else{
        //  printf("(%d, %d) not generated.\n", node->xpos, node->ypos-1);
            flag = 0;
            }
        }
        return;
    }

}

int point_in_pol(int vertcount, float *vertx, float *verty, int vertexx, int vertexy){
    double vertexx1;
    vertexx1 = vertexx;
    double vertexy1;
    vertexy1 = vertexy;
  int i ,j, c = 0;
  for (i = 0, j = vertcount-1; i < vertcount; j = i++) {
    if ( (((verty[i]>=vertexy1) && (verty[j]<=vertexy1) )  ||  ((verty[i]<=vertexy1)   && (verty[j]>=vertexy1) )) &&
     (vertexx1 < (vertx[j]-vertx[i]) * (vertexy1-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
   // if ((vertexx1 == ((vertx[j]-vertx[i]) * (vertexy1-verty[i]) / (verty[j]-verty[i])+ vertx[i])) )
//  return 1;
  }
  return c;
}

int main(){

    printf("\nFILE NUMBER 1\n");

    int x_start, y_start, x_end, y_end;
    clock_t begin, end;
    double time_spent;
    FILE *fp;
    fp = fopen("1.txt", "r");
    if (fp == NULL)
        printf("Error printing output. \n");
    else
    fscanf(fp, "(%d,%d)\n", &x_start, &y_start);
    fscanf(fp, "(%d,%d)\n", &x_end, &y_end);
    printf("Starting point is (%d, %d)\n", x_start, y_start);
    printf("Ending point is (%d, %d)\n", x_end, y_end);
    char temp3[255];
    char* source;
    int t;
    while(fgets(temp3, 255, fp)!= NULL){
        source = temp3;
        t = 0;
        printf("Polygon %d: ", polycount);
        while(sscanf(source, "(%f,%f)%*[^(]%n", &polygonx[polycount][vertcount], &polygony[polycount][vertcount], &t) == 2){
            printf("(%.2f,%.2f),",polygonx[polycount][vertcount], polygony[polycount][vertcount]);
            source+=t;
            vertcount++;
        }
        printf("\n");
        vertcounts[polycount] = vertcount;
        vertcount = 0;
        polycount++;
    }
    printf("Select a search algorithm:\n 1. BFS\n 2: DFS\n 3: A*");
    scanf("%d", &des);
    if (des == 1 || des == 2){
        begin = clock();
        Node* start = generate(x_start, y_start);
        enqueue(start);
        while(queuesize!=0){
        expand(front->Gnode, x_end, y_end);
        }
        end = clock();
        time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
        printf("Solution cost is: %d.\n", solutioncost);
        printf("Running time is %f.\n", time_spent);

        fclose(fp);
    }
}

using this as input for 1.txt:

(4,4)
(1,1)
(3,2),(2,2),(2,3),(3,3)

will yield no answer

Point *on* edge should be a rare case. Fortunately, checking if a point is on a line segment is a straightforward calculation (its only drawback is indeed the finite mathematical precision of computers). I would suggest creating a separate function for it and only check if point-in-poly returns `false`. You may want to add a global bounding box check before calling either of these functions (and if you have lots of polys to check, cache their bounds). — Jongware, Sep 27 '15 at 13:44

4pie0 · Accepted Answer · 2015-09-28T14:19:48.743

int
point_in_pol(int vertcount, float *vertx, float *verty,
int vertexx, int vertexy)
{
    double vertexx1;
    vertexx1 = vertexx;
    double vertexy1;
    vertexy1 = vertexy;

    int i ,j, c = 0;
    for (i = 0, j = vertcount-1; i < vertcount; j = i++)
    {
        if ( (((verty[i]>=vertexy1) && (verty[j]<=vertexy1) )
            || ((verty[i]<=vertexy1) && (verty[j]>=vertexy1))) // this checks
                                                               // whether y-coord
                                                               // i.e. vertexy1 is
                                                               // between edge's
                                                               // vertices
            && (vertexx1 < (vertx[j]-vertx[i]) 
                * (vertexy1-verty[i]) / (verty[j]-verty[i])
                                              + vertx[i]) )    // this checks
                                                               // whether x-coord
                                                               // i.e. vertexx1
                                                               // is to the left
                                                               // of the line

       c = !c;
  }
  return c;
}

The prototype of this algorithm is named pnpoly and it's explanation can be found here. One of it's limitations is that it can't handle points located exactly on the boundary, i.e. it cannot say when it is the case.

Point on a (Boundary) Edge

Pnpoly partitions the plane into points inside the polygon and points outside the polygon. Points that are on the boundary are classified as either inside or outside.

I think this should do the trick:

if (vertexx1 == (vertx[j]-vertx[i]) 
     * (vertexy1-verty[i]) / (verty[j]-verty[i])
                                   + vertx[i]) ) // this will check if vertexx1
                                                 // is equal to x-coordinate
                                                 // of what would have point on the
                                                 // edge with y=vertexy1 had
    return 1;

But you should be careful about floating point error. Computational roundoff error will make the result wrong. You can add epsilon comparison:

if (fabs(vertexx1 - (vertx[j]-vertx[i]) 
     * (vertexy1-verty[i]) / (verty[j]-verty[i])
                                   - vertx[i]) < EPSILON)
    return 1;

It did not do the trick. :( Take note that i changed the first inequality to at at least get the correct output on polygons with horizontal lines (like boxes). — gkhnavarro, Sep 27 '15 at 14:15
@tinky_winky on which types of arguments do you think these floating point errors will occur? I have only been taking all-integer inputs in this function, and even the simplest ones do not work. — gkhnavarro, Sep 27 '15 at 15:13
@bhounter both verty[j]-verty[i] and vertx[j]-vertx[i] are floating. Could you please add link to complete code with test example? — 4pie0, Sep 27 '15 at 15:21
Im using this function for a pathfinding algorithm just so the nodes wont expand. I use all polygon edges as integers converted to floats, so vertx and verty just contain values like 1.00 , 2.00 I will edit my post for the whole code — gkhnavarro, Sep 27 '15 at 15:33

score 0 · Answer 2 · answered Sep 27 '15 at 15:28

A polygon can be convex or concave. If you have a convex polygon, then each border separates the space in two halves. If you determine whether a point is on a "good place" from a border's perspective, and it is on the "wrong side", then the point is outside the polygon. If the point is in "good place" from all the borders' perspective of a convex polygon, then it is inside the convex polygon.

If you have a concave polygon, then you can define your concave polygon as a set of convex polygons. If the point is inside any of those convex polygons, then it is inside the concave polygon.

If you have more than two dimensions, then you must first check whether the point is inside the plane where the polygon is. If yes, then you must analyze it as described above. If it is not in the same plane, then it is not inside the polygon.

can you please explain further? Do i still cast rays to do this? — gkhnavarro, Sep 27 '15 at 16:55
Each border of a convex polygon is a subset of finite length of the line which holds the border. If that line is strictly between an arbitrary point of the polygon and the point you check, then the point is outside the polygon. If there is no border whose backbone line is strictly between an arbitrary point of the polygon and the point you are analyzing, then the point is inside the convex polygon. If your polygon is concave, you need to divide it to convex polygons, reducing your task to several similar, but simpler tasks. — Lajos Arpad, Sep 27 '15 at 18:44

Including edges in this pnpoly algorithm

2 Answers2