convex hull - Graham's scam Algorithm version-0.1

input：a finite set of two dimentional points P (the number of points n is more than 2)

output: the convex hull of P 

Pesudo:

1, sort P in Lexgrahical order

2, construct the upper hull UHLL:

        1' add two points on UHLL

        2' for i=3 to  n

                 add P[i] on UHLL

                 While( UHLL.Length > 2 && the last 3 points of UHLL do not make a right turn)

                            delete the middle point of the last 3 points of UHLL

3, construct the lower hull LHLL:

        1' add two points on LHLL

        2' for i=3 to  n

                 add P[i] on LHLL

                 While( LHLL.Length > 2 && the last 3 points of LHLL do not make a right turn)

                            delete the middle point of the last 3 points of LHLL

4, delete the first and last point of LHULL

5, merge UHLL and LHULL to  convex hull of P CHULL

6, output CHULL

------------------------------------------------------------------------------------------------------------

next is objectarx implementation( with visual studio 2010 + AutoCAD 2013 x64 +object 2013):

 

first, get the point set:

 1 /********************************************************************************/
 2 /* get points in selection set                                                  */
 3 /* reference: http://www.theswamp.org/index.php?topic=30971.msg365170#msg365170 */
 4 /********************************************************************************/
 5 void CModifyEnt::getPoints(AcGePoint3dArray &vec)
 6 {
 7     AcDbObjectIdArray ids;
 8     CResbufPtr pFilter (acutBuildList(RTDXF0,ACRX_T("POINT,LINE,LWPOLYLINE"),NULL));
 9     CSelectionSet sset;
10     sset.userSelect(pFilter);
11     sset.asArray(ids);
12     if (ids.length() < 3)
13     {
14         acutPrintf(L"there must be more than 3 entities!");
15         return;
16     }
17     for(int i = 0 ; i < ids.length(); i++)
18     {
19         AcDbEntityPointer pEnt(ids[i],AcDb::kForRead);
20         if(pEnt.openStatus() == eOk )
21         {
22             if(pEnt->isA() == AcDbPoint::desc())
23             {
24                 AcDbPoint *pPoint = AcDbPoint::cast(pEnt);
25                 if(pPoint)
26                 {
27                     vec.append(pPoint->position());
28                 }
29             }
30             else if(pEnt->isA() == AcDbLine::desc())
31             {
32                 AcDbLine *pLine = AcDbLine::cast(pEnt);
33                 if(pLine)
34                 {
35                     vec.append(pLine->startPoint());
36                     vec.append(pLine->endPoint());
37                 }
38             }
39             else if(pEnt->isA() == AcDbPolyline::desc())
40             {
41                 AcDbPolyline *pLine = AcDbPolyline::cast(pEnt);
42                 if(pLine)
43                 {
44                     for(unsigned int i = 0 ; i < pLine->numVerts() ; i++)
45                     {
46                         AcGePoint3d pt;
47                         pLine->getPointAt(i,pt);
48                         vec.append(pt);
49                     }
50                 }
51             }
52         }
53     }
54 }

 

second, sort points:

 

 1 void CCalculation::insertion_sort(AcGePoint3dArray &points)
 2 {
 3     int p;
 4     int k;
 5     AcGePoint3d current_element;
 6 
 7     int length = points.length();
 8     for ( p = 1; p < length; p++ )
 9     {
10         current_element = points.at(p);
11         k = p - 1;
12         while( k >= 0 && isLessOrEqual(current_element, points.at(k)))
13         {
14             points.at(k+1) = points.at(k);
15             k--;
16         }
17         points.at(k+1) = current_element;
18     }
19 }

 

in this method, there is a method called "isLessOrEqual" to compare two point in lexicographical order:

 1 bool CCalculation::isLessOrEqual(AcGePoint3d &p1, AcGePoint3d &p2)
 2 {
 3     if (p1.x == p2.x)
 4     {
 5         if (p1.y == p2.y)
 6         {
 7             return p1.z <= p2.z;
 8         }
 9         else
10         {
11             return p1.y < p2.y;
12         }
13     }
14     else
15     {
16         return p1.x < p2.x;
17     }
18 }

 

 

then construct the upper hull:

 

 1 void CCalculation::constructUpper()
 2 {
 3     this->m_upperHull.removeAll();
 4     this->m_upperHull.append(this->m_points.at(0));
 5     this->m_upperHull.append(this->m_points.at(1));
 6 
 7     int k;
 8     AcGeVector3d v1, v2;
 9     for (int i=2; i < this->m_points.length(); i++)
10     {
11         this->m_upperHull.append(this->m_points.at(i));
12 
13         while(this->m_upperHull.length() > 2)
14         {
15             k = this->m_upperHull.length() - 1;//the last point of current hull's index
16 
17             v1 = this->m_upperHull.at(k) - this->m_upperHull.at(k-2);
18             v2 = this->m_upperHull.at(k-1) - this->m_upperHull.at(k-2);
19 
20             if(CCalculation::isRightTurn(v1,v2))
21                 break;
22 
23             this->m_upperHull.removeAt(k-1);
24         }
25     }
26 }

 

 

 

then construct the lower hull:

 

 1 void CCalculation::constructLower()
 2 {
 3     this->m_lowerHull.removeAll();
 4     this->m_lowerHull.append(this->m_points.at(0));
 5     this->m_lowerHull.append(this->m_points.at(1));
 6 
 7     int k;
 8     AcGeVector3d v1, v2;
 9     for (int i=2; i < this->m_points.length(); i++)
10     {
11         this->m_lowerHull.append(this->m_points.at(i));
12 
13         while(this->m_lowerHull.length() > 2)
14         {
15             k = this->m_lowerHull.length() - 1;//the last point of current hull's index
16 
17             v1 = this->m_lowerHull.at(k) - this->m_lowerHull.at(k-2);
18             v2 = this->m_lowerHull.at(k-1) - this->m_lowerHull.at(k-2);
19 
20             if(CCalculation::isLeftTurn(v1,v2))
21                 break;
22 
23             this->m_lowerHull.removeAt(k-1);
24         }
25     }
26 }

 

in this method, there are three methods:

isRightTurn( if two vectors are in right turn order);

1 bool  CCalculation::isRightTurn(AcGeVector3d &p1, AcGeVector3d &p2)
2 {
3     return     crossProduction(p1, p2).z > 0;
4 }

 

isLeftTurn( if two vectors are in left turn order); 

1 bool  CCalculation::isLeftTurn(AcGeVector3d &p1, AcGeVector3d &p2)
2 {
3     return     crossProduction(p1, p2).z < 0;
4 }

 

crossProduction(calculate two vectors' cross production):

 

1 AcGeVector3d CCalculation::crossProduction(AcGeVector3d &p1, AcGeVector3d &p2)
2 {
3      return AcGeVector3d((p1.y*p2.z - p1.z*p2.y), (p1.z*p2.x - p1.x*p2.z), (p1.x*p2.y - p1.y*p2.x));
4 }

 

 

then construct the convex hull:

 

 1 AcGePoint3dArray CCalculation::getConvexHull()
 2 {
 3     this->initPoints();
 4     this->constructUpper();
 5     this->constructLower();
 6 
 7     this->m_convexHull.removeAll();
 8     for (int i=0; i < this->m_upperHull.length(); i++)
 9     {
10        this->m_convexHull.append(this->m_upperHull[i]);
11     }
12 
13     for (int i = this->m_lowerHull.length() - 2; i >=0 ; i--)
14     {
15         this->m_convexHull.append(this->m_lowerHull[i]);
16     }
17 
18     return this->m_convexHull;
19 }

 

 

the CCalculation.h:

 1 class CCalculation
 2 {
 3 public:
 4     CCalculation(void);
 5     ~CCalculation(void);
 6 
 7     static bool isLessOrEqual(AcGePoint3d &p1, AcGePoint3d &p2);
 8     static void insertion_sort(AcGePoint3dArray &points);
 9     static AcGeVector3d crossProduction(AcGeVector3d &p1, AcGeVector3d &p2);
10     static bool isRightTurn(AcGeVector3d &, AcGeVector3d &);
11     static bool isLeftTurn(AcGeVector3d &, AcGeVector3d &);
12     static void getPoints(AcGePoint3dArray &points);
13     void initPoints();
14     AcGePoint3dArray getUpperHull();
15     AcGePoint3dArray getLowerHull();
16     AcGePoint3dArray getConvexHull();
17 private:
18     AcGePoint3dArray m_points;//selected points
19     AcGePoint3dArray m_upperHull;//upper hull
20     AcGePoint3dArray m_lowerHull;//lower hull
21     AcGePoint3dArray m_convexHull;//convex hull
22 
23     void constructUpper();
24     void constructLower();
25     void constructConveHull();
26 };

 

the CCalculation.cpp:

  1 /********************************************************************************/
  2 /* get points in selection set                                                  */
  3 /* reference: http://www.theswamp.org/index.php?topic=30971.msg365170#msg365170 */
  4 /********************************************************************************/
  5 void CCalculation::getPoints(AcGePoint3dArray &points)
  6 {
  7     AcDbObjectIdArray ids;
  8     CResbufPtr pFilter (acutBuildList(RTDXF0,ACRX_T("POINT,LINE,LWPOLYLINE"),NULL));
  9     CSelectionSet sset;
 10     sset.userSelect(pFilter);
 11     sset.asArray(ids);
 12     if (ids.length() < 3)
 13     {
 14         acutPrintf(L"there must be more than 3 entities!");
 15         return;
 16     }
 17     for(int i = 0 ; i < ids.length(); i++)
 18     {
 19         AcDbEntityPointer pEnt(ids[i],AcDb::kForRead);
 20         if(pEnt.openStatus() == eOk )
 21         {
 22             if(pEnt->isA() == AcDbPoint::desc())
 23             {
 24                 AcDbPoint *pPoint = AcDbPoint::cast(pEnt);
 25                 if(pPoint)
 26                 {
 27                     points.append(pPoint->position());
 28                 }
 29             }
 30             else if(pEnt->isA() == AcDbLine::desc())
 31             {
 32                 AcDbLine *pLine = AcDbLine::cast(pEnt);
 33                 if(pLine)
 34                 {
 35                     points.append(pLine->startPoint());
 36                     points.append(pLine->endPoint());
 37                 }
 38             }
 39             else if(pEnt->isA() == AcDbPolyline::desc())
 40             {
 41                 AcDbPolyline *pLine = AcDbPolyline::cast(pEnt);
 42                 if(pLine)
 43                 {
 44                     for(unsigned int i = 0 ; i < pLine->numVerts() ; i++)
 45                     {
 46                         AcGePoint3d pt;
 47                         pLine->getPointAt(i,pt);
 48                         points.append(pt);
 49                     }
 50                 }
 51             }
 52         }
 53     }
 54 }
 55 
 56 void CCalculation::initPoints()
 57 {
 58     CCalculation::getPoints(this->m_points);
 59     CCalculation::insertion_sort(this->m_points);
 60 }
 61 
 62 bool CCalculation::isLessOrEqual(AcGePoint3d &p1, AcGePoint3d &p2)
 63 {
 64     if (p1.x == p2.x)
 65     {
 66         if (p1.y == p2.y)
 67         {
 68             return p1.z <= p2.z;
 69         }
 70         else
 71         {
 72             return p1.y < p2.y;
 73         }
 74     }
 75     else
 76     {
 77         return p1.x < p2.x;
 78     }
 79 }
 80 
 81 void CCalculation::insertion_sort(AcGePoint3dArray &points)
 82 {
 83     int p;
 84     int k;
 85     AcGePoint3d current_element;
 86 
 87     int length = points.length();
 88     for ( p = 1; p < length; p++ )
 89     {
 90         current_element = points.at(p);
 91         k = p - 1;
 92         while( k >= 0 && isLessOrEqual(current_element, points.at(k)))
 93         {
 94             points.at(k+1) = points.at(k);
 95             k--;
 96         }
 97         points.at(k+1) = current_element;
 98     }
 99 }
100 
101 AcGeVector3d CCalculation::crossProduction(AcGeVector3d &p1, AcGeVector3d &p2)
102 {
103      return AcGeVector3d((p1.y*p2.z - p1.z*p2.y), (p1.z*p2.x - p1.x*p2.z), (p1.x*p2.y - p1.y*p2.x));
104 }
105 
106 bool  CCalculation::isRightTurn(AcGeVector3d &p1, AcGeVector3d &p2)
107 {
108     return     crossProduction(p1, p2).z > 0;
109 }
110 
111 bool  CCalculation::isLeftTurn(AcGeVector3d &p1, AcGeVector3d &p2)
112 {
113     return     crossProduction(p1, p2).z < 0;
114 }
115 
116 void CCalculation::constructUpper()
117 {
118     this->m_upperHull.removeAll();
119     this->m_upperHull.append(this->m_points.at(0));
120     this->m_upperHull.append(this->m_points.at(1));
121 
122     int k;
123     AcGeVector3d v1, v2;
124     for (int i=2; i < this->m_points.length(); i++)
125     {
126         this->m_upperHull.append(this->m_points.at(i));
127 
128         while(this->m_upperHull.length() > 2)
129         {
130             k = this->m_upperHull.length() - 1;//the last point of current hull's index
131 
132             v1 = this->m_upperHull.at(k) - this->m_upperHull.at(k-2);
133             v2 = this->m_upperHull.at(k-1) - this->m_upperHull.at(k-2);
134 
135             if(CCalculation::isRightTurn(v1,v2))
136                 break;
137 
138             this->m_upperHull.removeAt(k-1);
139         }
140     }
141 }
142 
143 void CCalculation::constructLower()
144 {
145     this->m_lowerHull.removeAll();
146     this->m_lowerHull.append(this->m_points.at(0));
147     this->m_lowerHull.append(this->m_points.at(1));
148 
149     int k;
150     AcGeVector3d v1, v2;
151     for (int i=2; i < this->m_points.length(); i++)
152     {
153         this->m_lowerHull.append(this->m_points.at(i));
154 
155         while(this->m_lowerHull.length() > 2)
156         {
157             k = this->m_lowerHull.length() - 1;//the last point of current hull's index
158 
159             v1 = this->m_lowerHull.at(k) - this->m_lowerHull.at(k-2);
160             v2 = this->m_lowerHull.at(k-1) - this->m_lowerHull.at(k-2);
161 
162             if(CCalculation::isLeftTurn(v1,v2))
163                 break;
164 
165             this->m_lowerHull.removeAt(k-1);
166         }
167     }
168 }
169 
170 
171 
172 AcGePoint3dArray CCalculation::getUpperHull()
173 {
174     this->constructUpper();
175     return this->m_upperHull;
176 }
177 
178 AcGePoint3dArray CCalculation::getLowerHull()
179 {
180     this->constructLower();
181     return this->m_lowerHull;
182 }
183 
184 AcGePoint3dArray CCalculation::getConvexHull()
185 {
186     this->initPoints();
187     this->constructUpper();
188     this->constructLower();
189 
190     this->m_convexHull.removeAll();
191     for (int i=0; i < this->m_upperHull.length(); i++)
192     {
193        this->m_convexHull.append(this->m_upperHull[i]);
194     }
195 
196     for (int i = this->m_lowerHull.length() - 2; i >=0 ; i--)
197     {
198         this->m_convexHull.append(this->m_lowerHull[i]);
199     }
200 
201     return this->m_convexHull;
202 }

the getPoints method use two class(reference http://www.theswamp.org/index.php?topic=30971.msg365170#msg365170):

ResbufPtr.h:

 1 #pragma once
 2 #include "StdAfx.h"
 3 // Chuck Gabriel @ TheSwamp
 4 
 5 class CResbufPtr {
 6     resbuf* pHead; // Maintain a pointer to the head of the resource chain for de-allocation
 7     resbuf* pBuffer;
 8 public:
 9     CResbufPtr(resbuf* ptr) : pHead(ptr), pBuffer(ptr) {}
10     ~CResbufPtr() {
11         acutRelRb(pHead); // release the buffer
12         pHead = pBuffer = 0;
13     }
14     resbuf* operator->() { return pBuffer; }  // so you can do things like buf->resval.rstring
15     operator resbuf*() { return pBuffer; } // so you can pass the smart pointer to functions that expect a resbuf*
16     bool isNull() { return pBuffer == 0; } // null pointer check
17     void start() { pBuffer = pHead; } // in case you need to return to the head of the resource chain
18     void next() { pBuffer = pBuffer->rbnext; } // try to move the pointer to the next item in the resource chain
19 };

 

 SelectionSet.h

 1 #pragma once
 2 #include "arxHeaders.h"
 3 
 4 class CSelectionSet
 5 {
 6     ads_name ssname;
 7 public:
 8     CSelectionSet(void);
 9     ~CSelectionSet(void);
10     bool isNull();
11     Acad::PromptStatus userSelect(const resbuf *pFilter);
12     Acad::ErrorStatus asArray(AcDbObjectIdArray &setIds);
13 };

 SelectionSet.cpp

 1 #include "StdAfx.h"
 2 #include "SelectionSet.h"
 3 
 4 
 5 CSelectionSet::CSelectionSet(void)
 6 {
 7     ssname[0] = 0L;
 8     ssname[1] = 0L;
 9 }
10 
11 CSelectionSet::~CSelectionSet(void)
12 {
13     acedSSFree(ssname);
14 }
15 
16 
17 bool CSelectionSet::isNull()
18 {
19     return ssname[0] == 0L && ssname[1] == 0L;
20 }
21 
22 Acad::PromptStatus CSelectionSet::userSelect(const resbuf *pFilter)
23 {
24     acedSSFree(ssname);
25     return (Acad::PromptStatus) acedSSGet(NULL,NULL,NULL,pFilter,ssname);
26 }
27 
28 Acad::ErrorStatus CSelectionSet::asArray( AcDbObjectIdArray &setIds )
29 {
30     if(this->isNull())
31         return Acad::eNullPtr;
32     long sslength;
33     ads_name ent;
34     acedSSLength(ssname,&sslength);
35     for (long i = 0; i < sslength; i++)
36     {
37         AcDbObjectId oId;
38         acedSSName(ssname, i, ent);
39         if(acdbGetObjectId(oId, ent) == eOk)
40             setIds.append(oId);
41         else
42             return Acad::eNullObjectId;
43     }
44     return eOk;
45 }

the demo:

Done!

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

references:

1. http://www.theswamp.org/index.php?topic=30971.msg365170#msg365170

2. Computational Geometry - Algorithms and Applications(Mark de Berg · Otfried Cheong & Marc van Kreveld · Mark Overmars)

3. http://en.wikipedia.org/wiki/Graham_scan

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

 

posted @ 2014-09-19 14:42  大馋猫  阅读(804)  评论(0)  编辑  收藏  举报