A Matrix-Based 2-D Polygon Clipping Class

Quick shot of the demo application.

Introduction

This code is a two-dimensional polygon-clipping algorithm that determines precisely where a line intersects with a polygon border. This code works for both concave and convex polygons of completely arbitrary shape and is able to handle any line orientation. Furthermore, all arrays used to hold the line and polygon coordinates are allocated dynamically so the only limit on the number of polygon vertices and intersecting lines is the amount of available memory. Error checking is used to ensure that the bounds of the dynamic arrays are not overrun.

Background

This code was developed while working on dose calculation routines for cancer radiation treatment therapy. For these calculations, it was necessary to accurately determine the intersections of certain features, both anatomical and mechanical, with the radiation field. As a result, CPolygonClip was born. Due to its implementation in a clinical setting, this code is quick, accurate and robust. Some methods and functions have been left out due to time constraints, but will be added in future.

Theory or How It Works

This is the nasty mathematical portion of the article. You can, of course, use the class without reading this, but it is best to know what you are dealing with before you dive right in. Fortunately, most programmers will have the necessary mathematical prerequisites to understand the material.
In general, polygon clipping is a mathematical technique that evaluates whether 2 three-dimensional objects intersect one another. In the more limited context of this article, polygon-clipping is used to determine the points of intersection of a line segment with a polygon. To implement the CPolygonClip class, two assumptions have been made. First, the vertices of the polygon are arranged in a clock-wise fashion. Secondly, the coordinate system of the device context is assumed to be MM_TEXT, although it would be relatively simple to convert output from the class to any other system. To define a line segment, a starting point P₁(x₁,y₁) and an ending point P₂(x₂,y₂) are used, while the polygon is defined by a series of verticies V_i(X_i,Y_i), i = 1..n. (See figure.)

Polygon clipping illustration.

The coordinates of each point are considered the components of a vector whose z-component is unity. Now, let

S_i = determinant(V_i, P₁, P₂) or,

Matrix determinant.

S_i represents the equation of a line, the sign of which will be fixed for points lying on either side of the line. Therefore, by searching for a sign change in S_i between two consecutive vertices it is possible to determine if an intersection of line and polygon exists between the vertex pair. Finding a sign change will guarantee an intersection of the polygon with the line segment defined by P₁ and P₂. This calculation is repeated for all lines that may possibly intersect the polygon border. Once all values of S_i for all lines have been calculated, the total number of sign changes of S gives the total number of intersections of lines and polygon. Naturally, the total number of sign changes in S_i for a given line is indicative of the number of intersections that particular line makes with the polygon. S_i can also be used to determine if a line intersects precisely at a polygon vertex (collinear intersection). This situation is indicated by the condition S_i = 0. Having found the number of intersections, it is now possible to determine the coordinates of the intersection points. This is done by formulating a system of two linear equations in two unknowns. It is trivial to solve these and obtain the intersection point. We need two equations:

y = m_line*x + b_line

and

y = m_poly*x + b_poly

where y is the y coordinate, m_line the slope of the line, m_poly the slope of the polygon segment between the two vertices spanning the intersect, x the x-coordinate and b_line and b_poly the y-intercept of line and poly, respectively. These equations are easily solved by setting the y values equal - that is, setting one equation equal to another. This allows one to evaluate x, the x-intersection coordinate. Once x is known, either equation may be used to determine y. The coordinate point (x,y) gives the intersection point of line and polygon. Caution must be exercised for lines that are vertical or lines that have the same slope, since the equations above cannot be solved in these cases.

Using the code

Simply insert the .h and .cpp files into your current project and you are away. Call member functions as needed. The compiler will spew out some conversion warnings, but these are mainly due to the screen drawing calls (ie: double to int) and are not particularly significant when placing dots on a screen.

The class consists of a total of 7 functions and 4 adjustable parameters. As mentioned previously, every array used is allocated dynamically thus allowing any arbitrary polygon shape and number of lines to be used (subject to available memory). The CPolygonClip class is constructed using:

CPolygonClip(unsigned int nNumVerticies, unsigned int nNumLines);

Where the parameters passed are the number of vertices of the polygon and the number of lines for which to calculate the intersections. The number of vertices passed must be one greater than the actual number present in the polygon - the last vertex of the polygon must be a duplicate of the first vertex of the polygon.

unsigned int* CalcSiDeterm(double* LineX1, double* LineY1, 
			          double* LineX2, double* LineY2, 
			          double* PolyVertexX, double* PolyVertexY);

This function calculates the S_i determinants for all lines passed to it. LineX1 and LineY1 are pointers to arrays holding the x and y coordinate pairs of P₁, while LineX2 and LineY2 are pointers to arrays holding the x and y coordinate pairs of P₂. PolyVertexX and PolyVertexY are pointers to arrays holding the x and y coordinates of the polygon verticies. Again, it is important that these arrays have, as their last element, a vertex that is a duplicate of the first vertex of the polygon. CalcSiDeterm returns a pointer to an array that holds the number of intersections per line indexed by line number.

bool CPolygonClip::CalculateIntersections(double *LineX1, double *LineX2, 
double *LineY1, double *LineY2, 
double *PolyVertexX, double *PolyVertexY, 
unsigned int nMaxNumSignChanges)

CalculateIntersections() is fairly self-explanatory. It simply calculates the intersection coordinates of the points of intersection between the polygon and the line. It accepts the same parameters as CalcSiDeterm() with the exception of the last parameter, nMaxNumSignChanges which represents the maximum numbers of sign changes for all the given lines. This is returned by FindArrayMax() which is described below.

unsigned int CPolygonClip::FindArrayMax(unsigned int *pNumChanges)

FindArrayMax() simply accepts a pointer to an array and extracts the maximum value from that array. This output is passed as the parameter nMaxNumSignChanges to CalculateIntersections().

void CPolygonClip::Draw(CDC *pDC)

This function simply draws blue points at the location of the intersections. Finally,

double CPolygonClip::TurboDeterm(double Elem11, double Elem12, <BR>                                 double Elem13, double Elem21, double Elem22, <BR>                                 double Elem23)

is a helper function that calculates a quick determinant of a 3x3 matrix, accepting the first six elements of the matrix.

CheckValidPoint(CPoint *lpVertexPoints, int nCount, CPoint aPoint)

A private method to the class which prevents detection of points which to not lie on a finite line.

A typical usage of the class is as follows:

	//////////////////////////////////
	//    BEGIN POLYGON CLIPPING	
	//////////////////////////////////
	m_pPolygonClip = new CPolygonClip(NumberOfVerticies, NumberOfLines);
	ASSERT(m_pPolygonClip != NULL);

	unsigned int* pnSignChanges = 0;
	
	// Calculate the determinants
	pnSignChanges = m_pPolygonClip->CalcSiDeterm(pLineX1, pLineY1, <BR>                                  pLineX2, pLineY2, pPolyVertsX, pPolyVertsY);

       // Find maximum number of intersections
       unsigned int nNumChanges = m_pPolygonClip->FindArrayMax(pnSignChanges);

	// Calculate line/polygon intersections
	bool bSuccess = m_pPolygonClip->CalculateIntersections(pLineX1, <BR>                                       pLineX2, pLineY1, pLineY2, pPolyVertsX, <BR>                                       pPolyVertsY, nNumChanges);

	// should include some other error 
	// handling routine...try/catch, etc...
	// (this will do for now)
	ASSERT(bSuccess);		

	// Draw the intersection points
	m_pPolygonClip->Draw(pDC);
	//////////////////////////////////
	//    END POLYGON CLIPPING		
	//////////////////////////////////

Don't forget to call the destructor when finished with a particular instance of the class.

Points of Interest

A few things to note here. Unfortunately, due to the application this class was written for, it is necessary to call the member functions of the class in a certain order. There are flags present to ensure that calls are made in the correct sequence. This is more due to the architecture of the clinical software than anything else.
Secondly, the dynamic memory allocation was a bit of a chore - especially coming from a FORTRAN background where arrays are indexed from 1. In future, I would probably utilize a 2-D array class for this in order to make life easier.
As it stands, this is not the full implementation of the class. For certain reasons, I had to extract portions of the code in order to make it available. Therefore, some functions have not yet been implemented since they still need trimming. Included among these are functions to determine which segments of the line lie inside/outside of the polygon border and a function to dump some of the information to ASCII file.
The matrix method was used in favour of the vector algebraic method as I feel it is more elegant and matricies are far more suited to computer implementation than vector algebra.
At any rate, I think this provides a quick and easy to use implementation of some basic polygon clipping, routines for which I have not been able to find readily on the internet.

The demo app is fairly static, simply showing a polygon (in red) and some black lines. The blue dots indicate the points of intersection. However, the calculations are all performed at runtime (having placed all the working code in OnDraw()) so this is not simply some pre-drawn picture. Feel free to play around with the polygon verticies and line orientations to see that it works.

Update March 13, 2003

Discovered that for finite lines, points that did not actually lie on the line were being detected and intersection points were being calculated. To correct this a private member function was added that allows detection of points that lie solely on the line segment rather than on an infinite line.