Calculation of polygon border width

[Bronislaw]

Hi people,
does anyone know how to calculate the inner border coordinates of a polygon in cpp, if you want a fixed border width?

It should look like this:
https://i.stack.imgur.com/GUzIo.gif

I can't use a fixed distance from the outer edges along the anglebisectors, because the smaller the angle is, the longer the distance from the outer edge to the inner border edge should be for it to look like on the given link, otherwise (with fixed distance) it would look very distorted.

For example I've got a triangle with the following coordinates:

OuterCoordA=(170,20), OuterCoordB=(100,700), OuterCoordC=(172,550)

anglebisectorA = 1.6194, anglebisectorB = -1.29575, anglebisectorC = 3.36187;

Usually I would calculate the inner coordinate position like below, but the right distance depends on the angle size:

distance = ?

InnerCoordA(x,y) = OuterCoordA.x+distance*cos(anglebisectorA),OuterCoordA.y+distance*sin(anglebisectorA)

What I need is the right distance or the inner coordinate values. Is there a simple/short mathematic concept/formula to get those coordinates?
Many thx for any help.

[Duthomhas]

[edit]
Yeah, ignore this crap. See lastchance’s answer below. I don’t know what was wrong with me when I wrote this...
[/edit]

Yes, it is simple and short math, but in a computer it is a little more to program it.

I would recommend the following procedure as the cleanest:

For every point (with two intersecting segments), shown in your graphic as the black-rimmed green spots, you must find the normal. There are a number of simple ways to accomplish this, but geometrically, you just bisect the angle.

The new point will scale along this bisection ray. Here is a crude ASCII graphic:

           \ new point
            O············
  \    x  .`       ↑
   \    .`         d
    \ .`           ↓
     O···················
    outer edge point

You also need to know whether you are scaling in or out. Do this easily by tracking the side that has the most acute angles.

Finally, as to your question, you wish to know how far to make distance x so that distance d is inviolate.

For that you’ll need a little geometry. Do you see the right triangle there?

             O
       x   .`|
         .` _| d
       .`  | |
      O·······
    a  

You know angle a, the square angle, and length d. Use the arcsine to calculate x.

Finally, you only need to scale your point along the ray x units.


Alas, it is a bit of work to code this. Sorry.

[lastchance]

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cmath>
using namespace std;

const double SMALL = 1.0e-10;

struct Point{ double x, y; };

int main()
{
   //=====================================================//
   // Assumes:                                            //
   //    N+1 vertices:   p[0], p[1], ... , p[N-1], p[N]   //
   //    Closed polygon: p[0] = p[N]                      //
   //=====================================================//
   vector<Point> p = { { 0, 0 }, { 1, 0 }, { 2, 3 }, { 0, 1 }, { 0, 0 } };    // Outside vertices

   int N = p.size() - 1;
   vector<Point> q(N+1);               // This will hold the inside vertices
   double thickness = 0.1;             // Border thickness

   double a, b, A, B, d, cross;

   // Unit vector (a,b) along last edge
   a = p[N].x - p[N-1].x;
   b = p[N].y - p[N-1].y;
   d = sqrt( a * a + b * b );  
   a /= d;
   b /= d; 

   // Loop round the polygon, dealing with successive intersections of lines
   for ( int i = 0; i < N; i++ )
   {
      // Unit vector (A,B) along previous edge
      A = a;
      B = b;
      // Unit vector (a,b) along next edge
      a = p[i+1].x - p[i].x;
      b = p[i+1].y - p[i].y;
      d = sqrt( a * a + b * b );
      a /= d;
      b /= d;
      // New vertex
      cross = A * b - a * B;
      if ( abs( cross ) < SMALL )      // Degenerate cases: 0 or 180 degrees at vertex
      {
         q[i].x = p[i].x - thickness * b;
         q[i].y = p[i].y + thickness * a;
      }
      else                             // Usual case
      {
         q[i].x = p[i].x + thickness * ( a - A ) / cross;
         q[i].y = p[i].y + thickness * ( b - B ) / cross;
      }
   }

   // Close the inside polygon
   q[N] = q[0];


   // Write to file in order to plot by columns
   ofstream out( "output.txt" );
   #define SP << " " << fixed << setw( 12 ) <<
   for ( int i = 0; i <= N; i++ ) out SP p[i].x SP p[i].y SP q[i].x SP q[i].y << '\n';
   out.close();
}

https://imgur.com/a/YABvs

Output file (output.txt):

     0.000000     0.000000     0.100000     0.100000
     1.000000     0.000000     0.927924     0.100000
     2.000000     3.000000     1.771175     2.629754
     0.000000     1.000000     0.100000     0.958579
     0.000000     0.000000     0.100000     0.100000

Assumes polygon traversed ANTICLOCKWISE. If you traverse it CLOCKWISE (or if you make thickness negative) you will get an OUTSIDE border.


In case you were dying to ask ...

- Vector equation of a straight line is

x = xv + t ( a, b )

where ( a, b ) is a unit vector, xv is any point on that line, t is any real number.

- The inward normal vector is then ( -b, a ) for an anticlockwise 90 degree turn.

- The vector equation of a PARALLEL line, distance thick away is then

x = xv + thick( -b, a ) + t ( a, b )

- The vector of the other internal line (unit vector ( A, B ) ) meeting at that vertex is

x = xv + thick( -B, A ) + u ( A, B )

- Solve for t and/or u and hence get the intersection of these lines, giving the internal vertex. After careful algebra and rearrangement it produces the expressions given in the code.

- In the degenerate case (180 degrees at the "vertex") the internal point is just

x = xv + thick( -b, a )

To determine whether the vertices are taken anticlockwise, take points in order 0,1,2,...,N and find the vector area. If it points in the positive z direction then the vertices are anticlockwise; if in the negative direction then clockwise.


You don't need to find any angles anywhere.

[Duthomhas]

You know, I knew that stuff, but my brain got stupid. I’m so sorry. (And ashamed.) lastchance’s answer is so much better.

[lastchance]

Here is an extended version which calls a separate border() function allowing you to do multiple polygons easily.

It also uses the signed area of the polygon to detect whether the points are in anticlockwise or clockwise order and acts appropriately, so that you don't have to worry about that.

The requirement for the polygon to be closed (i.e. last point = first point) is retained, however. It is an irritating special case to deal with otherwise.

Output should go to screen as well now, for testing in C++ shell.

Graphic of output: https://imgur.com/a/YABvs

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cmath>
using namespace std;

const double SMALL = 1.0e-10;

// Structures/classes
struct Point{ double x, y; };

// Prototypes
void border( const vector<Point> &p, double thickness, vector<Point> &q );
double signedArea( const vector<Point> &p );
void output( ostream &strm, const vector<Point> &p, vector<Point> &q );


//======================================================================


int main()
{
   vector<Point> p = { { 0, 0 }, { 1, 0 }, { 2, 3 }, { 0, 1 }, { 0, 0 } };    // Outside vertices
   double thickness = 0.1;                                                    // Border thickness
   vector<Point> q;                                                           // Inside vertices (eventually)

   // Compute internal vertices
   border( p, thickness, q );

   // Output to screen and/or file
   ofstream out( "output.txt" );
   output( cout, p, q );
   output(  out, p, q );
   out.close();

   cout << "\nSigned area is " << signedArea( p ) << '\n';
}


//======================================================================


void border( const vector<Point> &p, double thickness, vector<Point> &q )
{
   //=====================================================//
   // Assumes:                                            //
   //    N+1 vertices:   p[0], p[1], ... , p[N-1], p[N]   //
   //    Closed polygon: p[0] = p[N]                      //
   //    No zero-length sides                             //
   // Returns (by reference, as a parameter):             //
   //    Internal poly:  q[0], q[1], ... , q[N-1], q[N]   //
   //=====================================================//
   int N = p.size() - 1;
   q.resize(N+1);              

   double a, b, A, B, d, cross;

   double displacement = thickness;
   if ( signedArea( p ) < 0 ) displacement = -displacement;     // Detects clockwise order

   // Unit vector (a,b) along last edge
   a = p[N].x - p[N-1].x;
   b = p[N].y - p[N-1].y;
   d = sqrt( a * a + b * b );  
   a /= d;
   b /= d; 

   // Loop round the polygon, dealing with successive intersections of lines
   for ( int i = 0; i < N; i++ )
   {
      // Unit vector (A,B) along previous edge
      A = a;
      B = b;
      // Unit vector (a,b) along next edge
      a = p[i+1].x - p[i].x;
      b = p[i+1].y - p[i].y;
      d = sqrt( a * a + b * b );
      a /= d;
      b /= d;
      // New vertex
      cross = A * b - a * B;
      if ( abs( cross ) < SMALL )      // Degenerate cases: 0 or 180 degrees at vertex
      {
         q[i].x = p[i].x - displacement * b;
         q[i].y = p[i].y + displacement * a;
      }
      else                             // Usual case
      {
         q[i].x = p[i].x + displacement * ( a - A ) / cross;
         q[i].y = p[i].y + displacement * ( b - B ) / cross;
      }
   }

   // Close the inside polygon
   q[N] = q[0];
}


//======================================================================


double signedArea( const vector<Point> &p )
{
   double A = 0;
   //========================================================//
   // Assumes:                                               //
   //    N+1 vertices:   p[0], p[1], ... , p[N-1], p[N]      //
   //    Closed polygon: p[0] = p[N]                         //
   // Returns:                                               //
   //    Signed area: +ve if anticlockwise, -ve if clockwise //
   //========================================================//
   int N = p.size() - 1;
   for ( int i = 0; i < N; i++ ) A += p[i].x * p[i+1].y - p[i+1].x * p[i].y;
   A *= 0.5;
   return A;
}


//======================================================================


void output( ostream &strm, const vector<Point> &p, vector<Point> &q )
{
   int N = p.size() - 1;

   // Formatting (adapt to your needs)
   strm.setf( ios::fixed );                 // fixed or scientific
   strm.precision( 6 );                     // Decimal digits after point
   #define SP << " " << setw( 12 ) <<       // Spacer and field width

   for ( int i = 0; i <= N; i++ ) strm SP p[i].x SP p[i].y SP q[i].x SP q[i].y << '\n';
}

//======================================================================