1

Lets say we have a Circle with Center(float_x,float_y) and Radius float_r. This is located on a grid or array like plane and we want to find every Cell(int i,int j) that gets intersected by the circle.

These Cells should be ordered in an array according to their angle to the x axis (x axis is horizontal positive to the right, y vertical positive in up direction).

The problem is that the center of the circle doesn't have to be exactly in the middle of a cell, I guess this denies the use of "point" symmetry of a cell.

I thought about pretending to have a cell symmetric circle and iterating through 90° and adding the offset to each point i am calculating but i am unsure if that would yield good results.

I would appreciate any Help or Ideas

Thanks

EDIT:

The following Code is for finding every point in the first quadrant(x+,y+) i haven't tried it yet but i am pretty sure it will work. Can be applied to the other quadrants too but then x/y iteration order and direction has to be changed. Since we start with xmax/pt_y the points will be ordered by their angle.

As soon as i have tested this i will mark this question as solved.

float pt_x,pt_y are the circle middle coordinates

float searchradius is the radius of the circle

float map_info.scale is the size of one cell in the grid

int maxx=round((pt_x+searchradius)/map_info.scale) getting max possible cell

            j=round(pt_y/map_info.scale);
            for (i=maxx; i >round(pt_x/map_info.scale); i--) {
                while(true)
                {
                    Point.x=i;
                    Point.y=j;
                    Points.push_back(Point);
                    if(sqrt(pow(i*map_info.scale,2)+pow((j+1)*map_info.scale,2))<=searchradius)
                    {
                        j++;
                    }
                    else break;
                }
            }
  • A question like this is a tough sell around here. You're better off making an attempt at solving this yourself and then if you aren't satisfied with the results of your attempt, asking questions about the attempt. Right now you have a math problem, but you haven't made it a programming problem yet. – user4581301 Dec 08 '20 at 21:00
  • try google circle rasterization with subpixel accuracy,...is the grid regular? rectangular or square? What is the size of the grid cell? you could try naive approach (parametric circle equation) and then add missing pixels (neighbors with hamming distance > 1). And yes symmetry is out of question... Another option is to compute all intersection with grid lines and then reorder by "atan2" or by dividing into 90 deg zones and order by significant axis... – Spektre Dec 09 '20 at 09:20
  • @Spektre thanks for your answer! Yes, the grid is regular. It is a rectangular grid with dynamic size(to be precise its a grid of type nav_msgs:::OccupancyGrid from the ROS environment). The grid cell Size is a variable of this grid type and won't allways be the same. Thanks for the name of my problem! I tried to google this before and wasn't able to find any results, i'll give it a try again. The Idea to search for grid line intersections seems to be very good, that will be my next approach! – Tristan9497 Dec 09 '20 at 14:23
  • @user4581301 thanks for letting my know i am quiet new here so i don't have a lot experience how to ask questions. – Tristan9497 Dec 09 '20 at 14:26
  • @Tristan9497 I finaly finished editing of my answer ... looks like I finally handled all the edge cases and quirks related to this... You need to check up to 4 last points in the list not just 2 (2 are enough for interior of the section but on joints 4 are needed) ... – Spektre Dec 10 '20 at 13:36

2 Answers2

0

So lets test the intersection with grid lines computed on floats (as I suggested in comments). I would start with dividing the circle into 4 sections ("quadrants") like this:

"quadrants"

In each section only one axis is major (has always bigger or equal difference to the other one) so we can iterate it directly by incrementing (without creating holes in curve). The minor axis is then computed using circles equation:

y = sqrt( r*r - (x-x0)*(x-x0)  )
x = sqrt( r*r - (y-y0)*(y-y0)  )

so simply loop the major axis in order by 4x (not nested) for loops and add the points. However whenever the minor axis computed is not directly on the grid we need to add 2 points where the second point has decremented major axis. This can cause duplicate points so to avoid duplicates we should also check if added point is not already added to the end of the list. At the end of last section we should also check against the start of the list (or remove last few duplicate points) I think there would be up to 4 points like this max.

When put all this together in C++/VCL it look like this:

//---------------------------------------------------------------------------
// global grid
float gx0=0.0;      // cell origin
float gy0=0.0;
float gxs=30.0;     // cell size
float gys=20.0;
//---------------------------------------------------------------------------
//  grid -> screen     |  screen -> grid
//  sx = gx0 + gx*gxs  |  gx = (sx-gx0)/gxs
//  sy = gy0 + gy*gys  |  gy = (sy-gy0)/gys
//---------------------------------------------------------------------------
void draw_cell(TCanvas *scr,int x,int y,DWORD col)  // screen, cell, color
    {
    x=floor(gx0+(float(x)*gxs));
    y=floor(gy0+(float(y)*gys));
    scr->Pen->Color=clDkGray;
    scr->Brush->Color=TColor(col);
    scr->Rectangle(x,y,x+gxs,y+gys);
    }
//---------------------------------------------------------------------------
void draw_grid(TCanvas *scr,int xs,int ys,DWORD col)    // screen, its resolution, color
    {
    int x0,y0,x1,y1,x,y;
    // visible grid
    x0=floor((float( 0)-gx0)/gxs);
    y0=floor((float( 0)-gy0)/gys);
    x1= ceil((float(xs)-gx0)/gxs);
    y1= ceil((float(ys)-gy0)/gys);
    for (y=y0;y<=y1;y++)
     for (x=x0;x<=x1;x++)
      draw_cell(scr,x,y,col);
    }
//---------------------------------------------------------------------------
void draw_circle(TCanvas *scr,float cx,float cy,float r,DWORD col)  // screen, circle, color
    {
    int i,ix,iy;
    int *pnt,pnts=0;            // output list of cells { x0,y0, x1,y1, ... }
    float x,y,x0,y0,x1,y1,xx,yy,rr=r*r;
    // ranges where x or y is major axis
    x0=floor(cx-(r*sqrt(0.5)));
    y0=floor(cy-(r*sqrt(0.5)));
    x1=ceil (cx+(r*sqrt(0.5)));
    y1=ceil (cy+(r*sqrt(0.5)));
    // prepare list
    pnt=new int[4*(x1-x0+1)];
    // this adds point (px,py) if the pnt[pi] ... pnt[pi-6] in list is not the same as (px,py)
    #define pnt_add(pi,px,py)                                         \
        {                                                             \
        int e=1;                                                      \
        if (pi>=0)                                                    \
            {                                                         \
            if ((pi+1<pnts)&&((pnt[pi+0]==px)&&(pnt[pi+1]==py))) e=0; \
            if ((pi-1<pnts)&&((pnt[pi-2]==px)&&(pnt[pi-1]==py))) e=0; \
            if ((pi-3<pnts)&&((pnt[pi-4]==px)&&(pnt[pi-3]==py))) e=0; \
            if ((pi-5<pnts)&&((pnt[pi-6]==px)&&(pnt[pi-5]==py))) e=0; \
            }                                                         \
        if (e){ pnt[pnts]=px; pnts++; pnt[pnts]=py; pnts++; }         \
        }
    // rasterize "quadrants" in order so no sorting is needed later
    for (x=x0;x<x1;x++)
        {
        xx=x-cx; xx*=xx; y=rr-xx; if (y<0.0) continue;
        y=sqrt(y); iy=float(cy-y); ix=x;
        if (y-floor(y)>1e-10) pnt_add(pnts-2,ix-1,iy);
                              pnt_add(pnts-4,ix  ,iy);
        }
    for (y=y0;y<y1;y++)
        {
        yy=y-cy; yy*=yy; x=rr-yy; if (x<0.0) continue;
        x=sqrt(x); ix=float(cx+x); iy=y;
        if (x-floor(x)>1e-10) pnt_add(pnts-2,ix,iy-1);
                              pnt_add(pnts-4,ix,iy  );
        }
    for (x=x1;x>x0;x--)
        {
        xx=x-cx; xx*=xx; y=rr-xx; if (y<0.0) continue;
        y=sqrt(y); iy=float(cy+y); ix=x;
                              pnt_add(pnts-2,ix  ,iy);
        if (y-floor(y)>1e-10) pnt_add(pnts-4,ix-1,iy);
        }
    for (y=y1;y>y0;y--)
        {
        yy=y-cy; yy*=yy; x=rr-yy; if (x<0.0) continue;
        x=sqrt(x); ix=float(cx-x); iy=y;
                              pnt_add(pnts-2,ix,iy  );
        if (x-floor(x)>1e-10) pnt_add(pnts-4,ix,iy-1);
        }
    // check duplicity between pnt start and end
    for (i=pnts-8;i+1<pnts;i+=2)
     if (pnt[0]==pnt[i+0])
      if (pnt[1]==pnt[i+1])
        {
        pnts=i;
        break;
        }

    // render the list (continuously change color to show points order)
    for (i=0;i+1<pnts;i+=2) draw_cell(scr,pnt[i+0],pnt[i+1],0x010000*(25+((230*i)/pnts)));

    // ---- you can ignore all the rest of the code its my debug rendering -----
    // debug numbers
    scr->Font->Color=TColor(0x0000FF00);
    scr->Font->Height=gys;
    scr->Brush->Style=bsClear;
    for (i=0;i+1<pnts;i+=2)
        {
        x=gx0+float(pnt[i+0])*gxs;
        y=gy0+float(pnt[i+1])*gys;
        scr->TextOutA(x,y,i>>1);
        }
    scr->Brush->Style=bsSolid;
    // debug circle overlay (ignore this)
    if (1)
        {
        int x,y,rx,ry,w=8;
        x =floor(gx0+(cx*gxs));
        y =floor(gy0+(cy*gys));
        rx=floor(r*gxs);
        ry=floor(r*gys);
        scr->Pen->Color=clYellow;
        scr->Brush->Style=bsClear;
        scr->Ellipse(x-rx,y-ry,x+rx,y+ry);
        scr->MoveTo(x-w,y);
        scr->LineTo(x+w,y);
        scr->MoveTo(x,y-w);
        scr->LineTo(x,y+w);
        scr->Brush->Style=bsSolid;
        }

    // free the list you should store/copy/export it somwhere
    delete[] pnt;
    #undef pnt_add
    }
//---------------------------------------------------------------------------

This is the output (including my debug rendering):

preview

The result is in pnt[pnts] array holding the x,y grid coordinates in order ...

This is far from optimized for example:

Spektre
  • 49,595
  • 11
  • 110
  • 380
  • Hey, thanks for your answer. In the mean time i spend some time to get this to work myself i will post an answer to that next to yours. I'll take a look at your solution and am happy to accept that as the correct one, even if both work. If you have the time i would enjoy to hear, what you are thinking of my solution since it seems you have a lot more experience. I'll post it asap – Tristan9497 Dec 11 '20 at 13:47
  • @Tristan9497 from what I see in your deleted answer its almost the same... The difference is that I also handle special cases on crossings where should be 2 poitns added per iteration of major axis and also I test up to 4 points of the already finished path for duplicates so the result is points in order and without duplicates – Spektre Dec 11 '20 at 14:33
  • yes, i realised that, and this should be fixed now but i am struggling with the visualization right now. But you are right it should almost be the same – Tristan9497 Dec 11 '20 at 14:39
  • Ok, i posted my answer as well and accepted yours. Thanks again for your help! – Tristan9497 Dec 11 '20 at 15:13
0

Ok, second try. I have tried to limit the use of more complex math to a minimum but i am sure this can be optimized.

First the code part with the iterations;

        double minx=round((pt_x-searchradius)/map_info.scale);
        double maxx=round((pt_x+searchradius)/map_info.scale);
        double miny=round((pt_y-searchradius)/map_info.scale);
        double maxy=round((pt_y+searchradius)/map_info.scale);


        //since circle wont ever be point symetrical we need to iterate through all quadrants
        //this is caused since we wont shift the circlecenter to the middle of the nearest field otherwise symmetry could be used
        Points.clear();
        //first quadrant
        int i,j;
        j=round(pt_y/map_info.scale);
        for (i=maxx; i >round(pt_x/map_info.scale); i--) {
            while(checkcellincircle(i,j,pt_x,pt_y,map_info,Points))
            {
                j++;
            }
            j--;
        }
        //second quadrant
        i=round(pt_x/map_info.scale);
        for (j=maxy; j >round(pt_y/map_info.scale); j--) {
            while(checkcellincircle(i,j,pt_x,pt_y,map_info,Points))
            {
                i--;
            }
            i++;
        }
        //third quadrant
        j=round(pt_y/map_info.scale);
        for (i=minx; i<round(pt_x/map_info.scale); i++) {
            while(checkcellincircle(i,j,pt_x,pt_y,map_info,Points))
            {
                j--;
            }
            j++;
        }
        //fourth quadrant
        i=round(pt_x/map_info.scale);
        for (j=miny; j <round(pt_y/map_info.scale); j++) {
            while(checkcellincircle(i,j,pt_x,pt_y,map_info,Points))
            {
                i++;
            }
            i--;
        }    

Second the part that handles weather a cell should belong to the circle or not:

bool checkcellincircle(int &i, int &j,double &pt_x,double &pt_y, map_inf &map,std::vector<geometry_msgs::Point32> &Points)
{
    //this function checks the shortest distance from pt_x/pt_y to the given cell
    //if the distance is <= the searchradius the point gets added to the Pointsvector and the function returns true
    //else the function returns false
    double cminx,cminy,dx,dy;
    geometry_msgs::Point32 Point;
    if(i*map.scale>pt_x)
        dx=abs(pt_x-i*map.scale-map.scale/2);
    else if(i*map.scale<pt_x)
        dx=abs(pt_x-i*map.scale+map.scale/2);
    else dx=abs(pt_x-i*map.scale);
    if(j*map.scale>pt_y)
        dy=abs(pt_y-j*map.scale-map.scale/2);
    else if(j*map.scale<pt_y)
        dy=abs(pt_y-j*map.scale+map.scale/2);
    else dy=abs(pt_y-j*map.scale);

    //getting coordinates of closest point on cell edge
    if((abs(dx)>=abs(dy))||dy==0)
    {
        cminx=map.scale;
        cminy=(abs(dy)/abs(dx))*map.scale;
    }

    else if((abs(dx)<abs(dy))||dx==0)
    {
        cminy=map.scale;
        cminx=(abs(dx)/abs(dy))*map.scale;
    }
    if(sqrt(pow(dx-cminx,2)+pow(dy-cminy,2))<=searchradius)
    {
        Point.x=i;
        Point.y=j;
        Points.push_back(Point);
        return true;

    }
    else return false;
}

Finaly a Picture of the Result 1m rad with 1m grid and 0.05m map resolution:

rasterized circle