Polygoner

Innsidetest

Vi har svært ofte behov for å finne ut om et punkt er inne i et generelt polygon. Det kan være at vi peker på noe på skjermen som skal identifiseres eller flyttes eller manipuleres eller vi skal dra opp en popupmeny assosiert med det aktuelle polygonformede objektet. Oftest er problemstillingen slik at vi kan forenkle situasjonen til til innsidetester i rektangler, men vi må kunne ta mer generelle situasjoner.

En vanlig måte å gjøre dette på er som følger:

Vi tenker oss et svært langt linjestykke gjennom det aktuelle punktet. Lang betyr at den helt sikkert har sine endepunkter utenfor polygonet.
Vi finner alle skjæringer mellom denne linja og kanter i polygonet,
Dersom antall skjæringer på den ene siden, f.eks. til venstre for punktet, er odde, så er punktet inne i polygonet ellers er det utenfor. Dette kan vi slå fast siden vi starter utenfor polygonet og går inn ved første skjæring, ut ved annen, inn ved tredje osv.

Vi kan f.eks. velge en linje parallell med en av koordinataksene.

inside1

En svakhet med algoritmen er at vi ved numeriske avrundinger eller sammenligninger kan miste skjæringspunkter dersom linja går gjennom hjørner i polygonet eller faller sammen med kanter i polygonet.

Funksjonen nedenfor benytter seg av Intersectionfunksjonen som er vist i modulen Linjer. Den tar ikke spesielle hensyn til hjørner eller parallelle kanter. Den benytter dessuten en generell skjæringsalgoritme og tar ikke hensyn til den effektivitetsgevinsten som ligger i at den ene linja i skjæringsalgoritmen er parallell med x-aksen. Legg merke til at vi samler ikke opp skjæringspunkter. Vi bare holder styr på om vi har et odde eller like antall skjæringer til venstre for det aktuelle punktet,


BOOL Inside(CPoint p,CArray<CPoint,CPoint&>& pol)
{
  // Returns TRUE if p is inside pol
  // pol can be convex or concave
  // result is not interpreteable if pol has crossing lines
  if(pol.GetSize()<3)
    return FALSE;
  // close the polygon temporarily
  pol.Add(pol[0]); // sentinel

  // Count intersections to the left of p
  // odd is hit, even is miss
  int pix;
  CPoint p1,p2,ps;
  BOOL Inside=FALSE;

  p1.x=-2*WORLDSIZE;  // smaller than the smallest
  p1.y=p.y;
  p2.x=2*WORLDSIZE;  // bigger than the biggest
  p2.y=p.y;
  for(pix=0;pix< pol.GetSize()-1;pix++)
    if(Intersection(p1,p2,pol[pix],pol[pix+1],ps))
      if (ps.x < p.x ) Inside=!Inside;

  pol.RemoveAt(pol.GetUpperBound()); // remove sentinel
  return Inside;
}

Vi ser at algoritmen benytter seg av en konstant, WORLDSIZE, for å bestemme et linjestykke som er større enn polygonets utstrekning. Et alternativ er at vi kjenner eller finner polygonets utstrekning, omgivende rektangel.

Ordne en sverm

Utgangspunktet er en "sverm" av punkter samlet i en punktarray. Rekkefølgen er tilfeldig og dersom vi trekker linjer mellom punktene vil de gå på kryss og tvers. Vi ønsker å ordne punktene slik at de beskriver et polygon som ikke har kryssende linjer. A er den opprinnelige svermen og B er et ordnet polygon.

For å få grep på dette må vi betrakte vinkler mellom linjer. Vi går fram slik:

Velg et punkt som anker, P0
Sorter de resterende punktene slik at Vinkel(P0,Pi) < Vinkel(Po,Pi+1) for alle 0<i<n, der n er antall punkter i svermen.

Med Vinkel mener vi følgende:

Siden vi bare skal sammenligne vinkler trenger vi ikke bruke trigonometriske funksjoner for å beregne nøyaktige vinkler. Vi kan være fornøyd dersom vi finner en funksjon som er monotont stigende i alle fire kvadranter, tilsvarende 0,360 grader. I Sedgewicks bok om algoritmer finner vi en slik funksjon som vi kan anvende:


float theta(CPoint p1,CPoint p2)
{
  // calculates a pseudo angle for the line p1->p2
  // according to Sedgewick
  float t;
  int dx=p2.x-p1.x; int ax=abs(dx);
  int dy=p2.y-p1.y; int ay=abs(dy);
  if((dx==0)&&(dy==0))
    t=0.0f;
  else
    t=(1.0f*dy)/(1.0f*(ax+ay));
  if(dx<0)
    t=2-t;
  else
    if(dy<0)
      t=4+t;
  return t*90.0f;
}

Algoritmen baserer seg på et uttrykk for vinklen som er den ene kateten dividert på summen av absoluttverdiene av katetene: dy/(ax+ay).

Med denne hjelpefunksjonen kan vi gå løs på å ordne punktene i svermen:


int Arrange_Polygon(CArray<CPoint,CPoint&>& p)
{
  // reorganize the points in p such
  // that they constitute a polygon with
  // non-crossing edges
  if(p.GetSize()<=3)
    return p.GetSize();
  int ix,minix,M;
  float minangle;
  CPoint tP;
  // find an anchor: the lowest y-value
  minix=0;
  for(ix=1;ix<p.GetSize();ix++)
    if(p[ix].y < p[minix].y)
      minix=ix;

  // put the anchor as point 0
  tP=p[0];p[0]=p[minix];p[minix]=tP;
  // make a simple selectsort for
  // the rest of the swarm based on
  // quasi angles with startpoint
  for(M=1;M<p.GetSize()-1;M++)
  {
    minix=M;
    minangle=theta(p[0],p[M]);
    for(ix=M+1;ix<p.GetSize();ix++)
    {
      if(minangle>theta(p[0],p[ix]))
      {
        minangle=theta(p[0],p[ix]);
        minix=ix;
      }
    }
    tP=p[minix];p[minix]=p[M];p[M]=tP;
  }
  return p.GetSize();
}

Lage omsluttende polygon

Igjen er utgangspunktet en "sverm" av punkter uten noen meningsfull rekkefølge. Nå ønsker vi å finne et konvekst polygon som omslutter svermen. A er den opprinnelige svermen, og B er det omgivende polygonet.

Igjen bruker vi hjelpefunksjonen Theta i en rutine som samler det omgivnde polygonet i de nødendige første elementene i polygonet. Funksjonen er i likhet med Theta insprert av Sedgewick.:


int Convex_Hull(CArray<CPoint,CPoint&>& p)
{
  // reorganize the points in p such
  // that the first ones are the convex hull
  // according to Sedgewick
  // returns the last index of the hull
  if(p.GetSize()<=3)
    return p.GetSize();
  int ix,minix,M;
  float minangle,v;
  CPoint tP;
  // find a place to start: the lowest y-value
  minix=0;
  for(ix=1;ix<p.GetSize();ix++)
    if(p[ix].y < p[minix].y)
      minix=ix;

  M=-1;    // upper index for hullnodes
  p.Add(p[minix]);// sentinel
  minangle=0.0f;
  do
  {
    M+=1;
    // swap elements at M and min
    tP=p[M];p[M]=p[minix];p[minix]=tP;
    minix=p.GetUpperBound();
    v=minangle; minangle=360.0f;
    for(ix=M+1;ix<p.GetSize();ix++)
    {
      float teta=theta(p[M],p[ix]);
      if((teta > v)&&(teta < minangle))
      {
        minix=ix;
        minangle=teta;
      }
    }
  } while(minix!=p.GetUpperBound());
  p.RemoveAt(p.GetUpperBound());  // remove the sentinel
  return M;  // the M first points are the hull
}

Beregne areal

Vi ønsker å lage en arealberegning som fungerer på et vilkårlig polygon. Vi stiller bare krav til at polygonet ikke har kryssende kanter.

Algoritmen er basert på å vekselvis legge til eller trekke fra områder under kantene i polygonet.

Vi betrakter arealet under hvert av kantene relativt til basislinja under polygonet. Bidraget fra en kant er summen av et rektangel og en trekant som antydet til venstre på figuren:

(p[ix+1].x-p[ix].x)*min(p[ix+1].y,p[ix].y)+
0.5*(p[ix+1].x-p[ix].x)*abs(p[ix+1].y-p[ix].y)
=
 (p[ix+1].x-p[ix].x)*(min(p[ix+1].y,p[ix].y)+
0.5* abs(p[ix+1].y-p[ix].y))

Dersom vi regner oss rundt for alle kantene vil uttrykket (p[ix+1].x-p[ix].x) endre fortegn alt ettersom vi skal legge til eller trekke fra. Vi trenger ikke vite på forhånd om punktene er ordnet med eller mot sola, vi kan bare ta absoluttvedien til slutt.


float Area(CArray<CPoint,CPoint&>& pol)
{
  // calculates the area of the polygon pol
  // the result of a polygon with crossing
  // edges is not meaningfull

  if(pol.GetSize()<3)
    return 0.0f;

  // assume y=0 is lower than any point in pol
  // close the polygon temporarily
  pol.Add(pol[0]);// sentinel
  double a=0.0;
  // add or subtract a rectangle and a triangle
  for(int ix=0;ix<pol.GetSize()-1;ix++)
    a+=(pol[ix+1].x-pol[ix].x)*
      (min(pol[ix+1].y,pol[ix].y)
      +0.5*fabs((pol[ix+1].y-pol[ix].y)));
  // remove the sentinal
  pol.RemoveAt(pol.GetUpperBound());
  return (float)fabs(a);

}

Demo

Du kan leke med disse algoritmene i canvaset nedenfor. Hvis du skrur på testdata skal de grønne punktene ligge innenfor og de røde utenfor.

Algoritmene er skrevet om til Javascript i canvaset over.

/*---------------------------------------
make random points on a canvas: makeNewSet()
Producing polygon with no linecrossing: arrangePolygon()
Producing conveks polygon: produceHull()
Testing inside in both cases: Inside()
Calculation area as precentage of complete frame: Area()
Sedgewicks algoritms
bs 2014
-----------------------------------------*/
// the device content in the canvas
var dc=null;
// dimensions of the area we use (preferable whole canvas)
// see init() at bottom of file;
var W=0;
var H=0;
// set up limits
var radius=8;
var maxcount=100;
// actual points generated
var count=0; 
// actual points describing the hull (convex polygon)
var hullcount=0;
var HullIsSet=false;
var Arranged=false;
// global in inside test
var ps=null;
// coordinates for polygon
var x;
var y;
// coordinates for hit tests
var hitx=null;
var hity=null;
// to hold a point, only used in insidetest
function Point(x,y){
    this.x=x;
    this.y=y;
}
////// set up canvas and make a dataset /////
function init(cid,w,h)
{
    var canvas = document.getElementById(cid);
    if (canvas.getContext){
        W=w;
        H=h;
        dc = canvas.getContext('2d');
        makeNewSet();
   } 
   else 
        alert('cannot do canvas');
}
//// paint what ever is to be painted /////
function paint()
{
   if(dc==null){
      console.log("no dc");
      return;
   }
   var lim;
   dc.fillStyle = "rgb(255,255,255)";
   dc.fillRect(0,0,W,H);
   // points
   for(var ix=0;ix < count;ix++){
      dc.beginPath();
      dc.fillStyle="rgb(0,0,0)";
      dc.arc(x[ix],y[ix],radius/2,0,2*Math.PI,false);
      dc.fill();
   }
   // hull only or all
   var lim=count;
   if(HullIsSet)
      lim=hullcount;
   dc.strokeStyle="rgb(100,100,100)";
   dc.moveTo(x[0],y[0])
   for(var ix=0;ix < lim;ix++)
      dc.lineTo(x[ix],y[ix]);
   // close it
   dc.lineTo(x[0],y[0]);
   dc.stroke();
   
   // mark inside/outside of testpoints
   if((HullIsSet ||Arranged) && (hitx !=null)&&(hitx.length>0) )
   {
        for(var ix=0;ix < hitx.length;ix++){// a point is set
            dc.beginPath();
            if(Inside(new Point(hitx[ix],hity[ix]),lim))
                dc.fillStyle="rgb(0,200,0)";
            else
                dc.fillStyle="rgb(200,0,0)";
            dc.arc(hitx[ix],hity[ix],radius/4,0,2*Math.PI,false);
            dc.fill();
        }
    }
}
  
////////  make a new set of points ////////////
function makeNewSet()
{
    // randomize the number
    count=getRandom(5,maxcount);
    x=new Array();
    y=new Array();
    for(var ix=0;ix < count;ix++)
    {
        // random integer coordinates
        x[ix]=getRandom(radius,W-2*radius);
        y[ix]=getRandom(radius,H-2*radius);
    }
    HullIsSet=false;
    Arranged=false;
    paint();
}
///////// make hitpoints for inside test ////
function makeHitPoints(n){
    if(hitx==null){
        hitx=new Array();
        hity=new Array();
        for(var ix=0;ix < n;ix++)
        {
            // random floats in coordinate range
            hitx[ix]=1.0*getRandom(radius,W-2*radius)+1/Math.random();
            hity[ix]=1.0*getRandom(radius,H-2*radius)+1/Math.random();
        }
    }
    else{
     hitx=null;
     hity=null;
    }
    paint();
}

///////// make random integer /////
function getRandom(min,max){
    return Math.floor((Math.random() * max) + min); 
}

/////////// helper for rough angle-equiv.  //////////////
function theta(x1,y1,x2,y2)
{
        // calculates a pseudo angle for the line p1->p2
        // according to Sedgewick
        var t;
        var dx=x2-x1; var ax=Math.abs(dx);
        var dy=y2-y1; var ay=Math.abs(dy);
        if((dx==0)&&(dy==0))
            t=0.0;
        else
            t=(1.0*dy)/(1.0*(ax+ay));
        if(dx<0)
            t=2.0-t;
        else 
            if(dy<0)
                t=4.0+t;
        return t*90.0;
}
//////  swap coordinates //////////////////
function swap(ix1,ix2){
    var tmpx=x[ix1];  var tmpy=y[ix1]; 
    x[ix1]=x[ix2]; y[ix1]=y[ix2];
    x[ix2]=tmpx; y[ix2]=tmpy;
}
//////// make a convex hull ///////////////////
function makeHull()
{
    // reorganize the points in p such
    // that the first ones are the convex hull
    // according to Sedgewick
    if(count<=3)
        return;
    //int ix,minix,M,tmpx,tmpy;
    //float minangle,v;
    // find a place to start: the lowest y-value
    // which is necessarily a part of the convex hull border
    var minix=0;
    for(var ix=1;ix<count;ix++)
        if(y[ix] < y[minix])
            minix=ix;
    // close the curve
    // add sentinel
    x[count]=x[minix];
    y[count]=y[minix];
    var minangle=0.0;
    var M=-1;        // upper index for hullnodes
    do
    {
        M+=1;
        // swap elements at M and minix
        swap(M,minix);
        minix=count;
        var v=minangle; 
        // start hi
        minangle=360.1;
        for(var ix=M+1;ix < count+1;ix++)
        {
            var teta=theta(x[M],y[M],x[ix],y[ix]);
            if((teta > v) && (teta < minangle))
            {
                minix=ix;
                minangle=teta;
             }
         }
    } while(minix != count);

    hullcount=M+1;
    HullIsSet=true;
    Arranged=false;
    paint();
}
//////////// arrange polygon without crossings ////////
function arrangePolygon()
{
    // reorganize the points in x,y such
    // that they constitute a polygon with 
    // non-crossing edges
    if(count<=3)
        return;
    //int ix,minix,M,tmpx,tmpy;
    var minangle;
    // find an anchor: the lowest y-value
    var minix=0;
    for(ix=1;ix<count;ix++)
        if(y[ix] < y[minix])
            minix=ix;
    // put the anchor as point 0
    swap(0,minix);
   
    // make a simple selectsort for 
    // the rest of the swarm based on
    // quasi angles with startpoint
    for(M=1;M<count-1;M++)
    {
        minix=M;
        minangle=theta(x[0],y[0],x[M],y[M]);
        for(ix=M+1;ix<count;ix++)
        {
            if(minangle > theta(x[0],y[0],x[ix],y[ix]))
            {
                minangle=theta(x[0],y[0],x[ix],y[ix]);
                minix=ix;
            }
        }
        swap(minix,M);
    }
    HullIsSet=false;
    Arranged=true;
    paint();
}
/////// calculate area ///////////////////////
function Area()
{
    // calculates the area of the polygon as 
    // percentage of available space
    if(count<3)
        return 0.0;
    // assume y=0 is lower than any point in pol
    // close the polygon temporarily
    // x[count]=x[0];y[count]=y[0];// sentinel
    var a=0.0;
    // add or subtract a rectangle and a triangle
    var lim;
    if(!HullIsSet)
        lim=count;
    else
        lim=hullcount;
    for(var ix=0;ix<lim;ix++)
    {
        // close the path
        if(ix==lim-1)
            a+=(x[0]-x[ix])*
                (Math.min(y[0],y[ix])
                +0.5*Math.abs((y[0]-y[ix])));
        else
            a+=(x[ix+1]-x[ix])*
                (Math.min(y[ix+1],y[ix])
                +0.5*Math.abs((y[ix+1]-y[ix])));
    }
    return Math.abs(a)*100.0/(1.0*(W*H));
}

//////// intersection between lines ///////////////
function Intersection(p1,p2,p3,p4,ps)
{
  // finds the intersecton between lines p1-p2 and p3-p4
  // return TRUE if intersection, FALSE else
  // result in ps, if intersection
  // use parametric equations for the lines
    
  var dx1=1.0*(p2.x-p1.x);
  var dx2=1.0*(p4.x-p3.x);
  var dy1=1.0*(p2.y-p1.y);
  var dy2=1.0*(p4.y-p3.y);
  var n=1.0*(dx2*dy1-dy2*dx1);
  if(n==0) return false;
  var s=1.0*(dx1*(p3.y-p1.y)-dy1*(p3.x-p1.x))/(1.0*n);
  if ((s<0.0)||(s>1.0)) return false;
  var t=1.0*(dy2*(p1.x-p3.x)-dx2*(p1.y-p3.y))/(1.0*n);
  if((t<0.0)||(t>1.0)) return false;
  ps.x=(dx1*t+p1.x);
  ps.y=(dy1*t+p1.y);
  return true;
}
////// testing inside point in polygon //////////////
function Inside(p, npnts)
{
  // Have we hit the polygon described int arrays x,y (npnts) ?
  if(npnts<3)
     return false;
  // Count intersections to the left of p
  // odd is hit, even is miss
  var p1=new Point(-100000,p.y);
  var p2=new Point(100000,p.y);
  var isInside=false;
  for(var pix=0;pix< npnts;pix++)
  {
     var p3= new Point(x[pix],y[pix]);
     var p4= new Point(x[pix+1],y[pix+1]);
        if(pix==npnts-1)
            p4=new Point(x[0],y[0]);
        
        ps= new Point(0,0);
     if(Intersection(p1,p2,p3,p4,ps))
        if (ps.x < p.x ) 
                isInside=!isInside;
  }
  return isInside;
}
/////  where it all starts //////
init('poly',400,400);

Eksempelprogram (C++/Visual Studio): https://svnit.hiof.no/svn/psource/CppGrafikk/planalg

Referanser

AlgorithmsRobert Sedgewick1984Addison Wesley0-201-06672-6