Click here to Skip to main content
15,904,652 members
Articles / Desktop Programming / MFC

Classes for computational geometry

Rate me:
Please Sign up or sign in to vote.
5.00/5 (25 votes)
26 Dec 2001CPOL1 min read 268.6K   3.8K   117   43
Some classes and utility functions for general computational geometry


This article presents two classes and a set of utility functions for computational geometry. C3Point is a 3D counterpart to CPoint and CPolygon encapsulates a set of C3Point's and provides general polygon handling functions. The classes have been mildly optimised for speed. The classes were originally written for use in discretising 2D surfaces into element networks and for calculating properties of the resultant elements. Care must be taken when using some of the functions such as the curvature and area functions to ensure that the results returned by the functions are consistent with your needs and definitions.

The classes make use of a typedef REAL that is either double or float depending on whether USING_DOUBLE or USING_FLOAT has been defined in geometry.h. Obviously using template classes would have been neater, but these classes were developed to get a job done, not to be the epitome of structured programming. A number of conversion functions have been provided:

D2Real(x) (x)             // double => REAL
F2Real(x) (x)             // float => REAL
Real2D(x) (x)             // REAL => double
Real2F(x) ((float)(x))    // REAL => float
Int2Real(x) ((double)(x)) // int => REAL
Real2Int(x) ((int)(x))    // REAL => int
Real2Int(double d0)       // REAL => int (faster than a (cast)).

All the classes and utility functions are provided 'as-is'. I've been meaning to write this class up for a long time and figured it was best to at least post something than nothing at all.


C3Point is a 3D counterpart to CPoint. It contains 3 data members x,y and z and a set of functions for calculating properties, scaling, translating and for arithmetic operations.

class C3Point {
// Attributes
    REAL x,y,z;

    C3Point() {}                          // constructor
    C3Point(double x, double y, double z) // constructor
    REAL Length2()                        // length squared
    REAL Length()                         // length
    void Scale(REAL factor)               // scale by a factor
    void Normalise();                     // convert to a unit length
    void operator=(C3Point& P)            // assign
    C3Point operator-(C3Point P)          // subtract
    C3Point operator-()                   // unary -
    C3Point operator+(C3Point P)          // add
    C3Point operator+=(C3Point P)         // add +=
    C3Point operator-=(C3Point P)         // subtract -=
    REAL operator*(C3Point P)             // vector dot product
    C3Point operator*(REAL f)             // scalar product
    C3Point operator/(REAL f)             // scalar div
    C3Point operator*=(REAL f)            // scalar mult *=
    C3Point operator/=(REAL f)            // scalar div /=
    C3Point operator^(C3Point P)          // cross product
    BOOL operator==(C3Point& P);          // is equal to?
    BOOL operator!=(C3Point& P)           // is not equal to?

#define VECTOR C3Point


CPolygon encapsulates a set of C3Point's and provides general polygon handling functions.

CPolygon(int);                     // Construct with a preallocated number of points

BOOL Closed();                     // Is the polygon closed?
int GetSize()                      // Number of points

// is vertex 'index' between start,end inclusive?
BOOL InSpan(int start, int end, int index);       
// is vertex 'index' between start,end exclusive?
BOOL InSpanProper(int start, int end, int index); 

BOOL PointIn(C3Point P);           // Is point inside polygon?
BOOL SetSize(int);                 // Change size of polygon
void RemoveAll()                   // Empty polygon of all points
BOOL Trim(int, int);               // Trims polygon down so that points before 
                                   //   "Start" and after "End" are removed.    
                                   //   Start and End must be in the range 0..GetSize()-1
BOOL Close();                      // Make polygon closed
BOOL Add(C3Point);                 // Add point to polygon
BOOL SetAt(int nPos, C3Point& p);  // set vertex nPos as point p
void Delete(int);                  // Delete a vertex
BOOL InsertAt(int nPosition, C3Point P); // insert point P at pos nPosition (0..N-1)

void FreeExtra();                  // Free extra memory left over after deletes

int PointSeparation(int Point1, int Point2); // Distance (in pts) between 2 points
void Rationalise(int nAngle);      // Combines adjacent line segments if the angle between 
                                   //   them is less than nAngle (degrees).
REAL SegmentLength(int,int);       // Length of a segment of the polygon
C3Point GetClosestPoint(C3Point p, int *nIndex = NULL);
REAL Area();                       // returns polygon area
C3Point Centroid();                // Calculate centroid of polygon
BOOL GetAttributes(REAL *pArea, 
                   C3Point *pCentroid, 
                   C3Point *pNorm,
                   REAL *pSlope, 
                   REAL *pAspect);
BOOL Diagonal(int i, int j);             // Returns TRUE iff (v_i, v_j) is a proper 
                                         //   internal or external diagonal of this polygon
virtual void Translate(VECTOR v);        // Translate polygon
BOOL Intersected(C3Point& p1, C3Point& p2);     // Does p1-p2 intersect polygon?
BOOL IntersectedProp(C3Point& p1, C3Point& p2); // Does p1-p2 intersect polygon properly?
BOOL Triangulate(CPolygon*);            // Triangulate: Ear clip triangulation
BOOL CPTriangulate(CPolygon*, C3Point); // Central point triangulation
BOOL DelauneyTri(CPolygon*);            // Triangulate: Delauney triangulation
// Load polygon from X-Y or X-Y-Z data file
BOOL LoadXY(LPCTSTR Filename, REAL Zdefault = D2Real(0.0));
BOOL LoadXY(FILE* fp, REAL Zdefault = D2Real(0.0));
BOOL LoadXYZ(LPCTSTR Filename, BOOL bWarn = TRUE);

// Save file either as:
//    Num Points, elevation, x-y pairs...,
// or
//    x-y-z triplets...
BOOL Save(LPCTSTR Filename, BOOL bAsPoints = FALSE, BOOL bWarn = TRUE);

void NaturalSpline(double*& b, double*& c, double*& d); // Natural cubic spline
REAL Curvature(int i);                                  // Curvature at vertex i
REAL Curvature(int nIndex, int nSampleSize);            // Avg curvature at i over 
                                                        // a number of points

C3Point& operator[](int index);
C3Point& Point(int index);
void operator=(CPolygon& P);

General Functions

These functions provide general routines for vectors (C3Points) and polygons.

inline REAL Dot(C3Point V1, C3Point V2)       // dot product
inline C3Point Cross(C3Point p1, C3Point p2)  // cross product
C3Point GetClosestPoint2D(C3Point& start, C3Point& end, C3Point& P);
REAL   Angle(C3Point, C3Point, C3Point);    // Angle between 2 vectors formed from 3 points (deg)
REAL   Angle(VECTOR v, VECTOR u);           // Angle between 2 vectors (degrees)
REAL   TriArea2(C3Point, C3Point, C3Point); // Area^2 between 2 vectors formed from 3 points 
REAL   TriArea2(VECTOR u, VECTOR v);        // Area^2 between 2 vectors
BOOL   IntersectProp(C3Point a, C3Point b,  // Returns true iff ab properly intersects cd: 
                     C3Point c, C3Point d)  //    they share a point interior to both segments.
                                            //    The properness of the intersection is ensured 
                                            //    by using strict leftness.

BOOL   Intersect(C3Point a, C3Point b,             // Returns true iff segments ab and cd
                 C3Point c, C3Point d);            // intersect, properly or improperly.
BOOL   Left(C3Point a, C3Point b, C3Point c);      // Returns true iff c is strictly to the left 
                                                   //   of the directed line through a to b.
BOOL   LeftOn(C3Point a, C3Point b, C3Point c);    // Same as Left, but c may be on the line ab.
BOOL   Colinear(C3Point a, C3Point b, C3Point c);  // Returns TRUE if a,b,c are colinear
BOOL   Between(C3Point a, C3Point b, C3Point c);   // Returns TRUE iff (a,b,c) are collinear and 
                                                   //   point c lies on the closed segement ab.
VECTOR Normal(C3Point p1, C3Point p2, C3Point p3); // Computes the normal (NOT unit normal) of 
                                                   //   a triangle, with points in Counter 
                                                   //   Clockwise direction.
VECTOR Scale(REAL factor, VECTOR v);               // Scales a vector by a factor.


The algorithms used are based in part from the book Computational Geometry in C by Joseph O'Rourke.


This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Written By
Founder CodeProject
Canada Canada
Chris Maunder is the co-founder of CodeProject and, and has been a prominent figure in the software development community for nearly 30 years. Hailing from Australia, Chris has a background in Mathematics, Astrophysics, Environmental Engineering and Defence Research. His programming endeavours span everything from FORTRAN on Super Computers, C++/MFC on Windows, through to to high-load .NET web applications and Python AI applications on everything from macOS to a Raspberry Pi. Chris is a full-stack developer who is as comfortable with SQL as he is with CSS.

In the late 1990s, he and his business partner David Cunningham recognized the need for a platform that would facilitate knowledge-sharing among developers, leading to the establishment of in 1999. Chris's expertise in programming and his passion for fostering a collaborative environment have played a pivotal role in the success of Over the years, the website has grown into a vibrant community where programmers worldwide can connect, exchange ideas, and find solutions to coding challenges. Chris is a prolific contributor to the developer community through his articles and tutorials, and his latest passion project, CodeProject.AI.

In addition to his work with, Chris co-founded ContentLab and DeveloperMedia, two projects focussed on helping companies make their Software Projects a success. Chris's roles included Product Development, Content Creation, Client Satisfaction and Systems Automation.

Comments and Discussions

QuestionSlow Delauney Pin
YvesDaoust10-Jul-17 4:47
YvesDaoust10-Jul-17 4:47 
AnswerRe: Slow Delauney Pin
Chris Maunder10-Jul-17 4:58
cofounderChris Maunder10-Jul-17 4:58 
GeneralLeast Squares algorithm or Linear Regression algorithm Pin
andrewtruckle10-Mar-08 4:32
andrewtruckle10-Mar-08 4:32 
GeneralVS2005 Pin
andrewtruckle28-Feb-08 23:12
andrewtruckle28-Feb-08 23:12 
QuestionCPolygon::Centroid problem Pin
andrewtruckle9-Aug-07 22:07
andrewtruckle9-Aug-07 22:07 
AnswerRe: CPolygon::Centroid problem Pin
Chris Maunder10-Aug-07 2:58
cofounderChris Maunder10-Aug-07 2:58 
GeneralRe: CPolygon::Centroid problem Pin
andrewtruckle12-Aug-07 23:39
andrewtruckle12-Aug-07 23:39 
GeneralRe: CPolygon::Centroid problem Pin
Chris Maunder13-Aug-07 4:34
cofounderChris Maunder13-Aug-07 4:34 
GeneralRe: CPolygon::Centroid problem Pin
andrewtruckle13-Aug-07 8:28
andrewtruckle13-Aug-07 8:28 
GeneralRe: CPolygon::Centroid problem Pin
Andrew Phillips29-Sep-08 3:00
Andrew Phillips29-Sep-08 3:00 
GeneralRe: CPolygon::Centroid problem Pin
andrewtruckle29-Sep-08 3:03
andrewtruckle29-Sep-08 3:03 
GeneralRe: CPolygon::Centroid problem Pin
Chris Maunder29-Sep-08 3:56
cofounderChris Maunder29-Sep-08 3:56 
GeneralRe: CPolygon::Centroid problem Pin
andrewtruckle29-Sep-08 5:04
andrewtruckle29-Sep-08 5:04 
GeneralThanks so much! Pin
FredericSivignonPro29-Jan-07 22:59
FredericSivignonPro29-Jan-07 22:59 
GeneralNatural Spline Pin
Zodiacon14-Jul-05 5:35
Zodiacon14-Jul-05 5:35 

The NaturalSpline function seems to give odd results sometimes.
Check out the following:
CPolygon p;
p.Add(C3Point(1875, 930, 0));
p.Add(C3Point(1868, 618, 0));
p.Add(C3Point(1990, 465, 0));

double *b, *c, *d;
p.NaturalSpline(b, c, d);

for(double x1 = 1875; x1 >= 1868; x1 -= 1) {
double x = x1 - 1875;
double y = x * x * x * d[0] + x * x * c[0] + b[0] * x + 930;
cout << "x=" << x1 << " y=" << y << "\n";

cout << "\n";

for(double x1 = 1868; x1 <= 1991; x1 += 4) {
double x = x1 - 1868;
double y = x * x * x * d[1] + x * x * c[1] + b[1] * x + 618;
cout << "x=" << x1 << " y=" << y << "\n";

delete []b; delete []c; delete []d;

The results for the second line segment seem way off:
Here's the output:
x=1875 y=930
x=1874 y=886.795
x=1873 y=843.419
x=1872 y=799.701
x=1871 y=755.471
x=1870 y=710.558
x=1869 y=664.792
x=1868 y=618

x=1868 y=618
x=1872 y=797.984
x=1876 y=959.468
x=1880 y=1103.08
x=1884 y=1229.44
x=1888 y=1339.19
x=1892 y=1432.95
x=1896 y=1511.34
x=1900 y=1574.99
x=1904 y=1624.53
x=1908 y=1660.59
x=1912 y=1683.8
x=1916 y=1694.77
x=1920 y=1694.15
x=1924 y=1682.55
x=1928 y=1660.6
x=1932 y=1628.93
x=1936 y=1588.17
x=1940 y=1538.94
x=1944 y=1481.87
x=1948 y=1417.59
x=1952 y=1346.73
x=1956 y=1269.9
x=1960 y=1187.75
x=1964 y=1100.9
x=1968 y=1009.96
x=1972 y=915.584
x=1976 y=818.38
x=1980 y=718.982
x=1984 y=618.017
x=1988 y=516.11

Although the final number converges to the third point, the intermmediate points create a very nasty looking curve.
Can you help?
GeneralIntersected only for 2D... Pin
Karl Runmo10-Dec-04 3:38
Karl Runmo10-Dec-04 3:38 
Questioncoordinate of intersection of line segments? Pin
Richard Males30-Aug-04 11:30
Richard Males30-Aug-04 11:30 
GeneralSubtract CPolygon from CPolygon Pin
anderslundsgard19-Jul-04 23:40
anderslundsgard19-Jul-04 23:40 
QuestionWhere is the colinear method? Pin
hgjgj121z20-Feb-04 12:46
hgjgj121z20-Feb-04 12:46 
GeneralGreat Job! Pin
Member 815358123-Jan-04 15:58
Member 815358123-Jan-04 15:58 
QuestionNon-MFC? Pin
TCMaster26-Aug-03 7:53
TCMaster26-Aug-03 7:53 
AnswerRe: Non-MFC? Pin
John M. Drescher26-Aug-03 8:59
John M. Drescher26-Aug-03 8:59 
GeneralRe: Non-MFC? Pin
TCMaster26-Aug-03 10:07
TCMaster26-Aug-03 10:07 
GeneralRe: Non-MFC? Pin
Antony M Kancidrowski28-Aug-03 2:01
Antony M Kancidrowski28-Aug-03 2:01 
GeneralRe: Non-MFC? Pin
Abin7-Sep-03 2:38
Abin7-Sep-03 2:38 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.