首页 > 代码库 > convex hull - Graham's scam Algorithm version-0.1
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 else10 {11 return p1.y < p2.y;12 }13 }14 else15 {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 index16 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 index16 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):
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 points19 AcGePoint3dArray m_upperHull;//upper hull20 AcGePoint3dArray m_lowerHull;//lower hull21 AcGePoint3dArray m_convexHull;//convex hull22 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 index131 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 index158 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 buffer12 pHead = pBuffer = 0;13 }14 resbuf* operator->() { return pBuffer; } // so you can do things like buf->resval.rstring15 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 check17 void start() { pBuffer = pHead; } // in case you need to return to the head of the resource chain18 void next() { pBuffer = pBuffer->rbnext; } // try to move the pointer to the next item in the resource chain19 };
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 else42 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
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
convex hull - Graham's scam Algorithm version-0.1
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。