设为首页收藏本站

中国膜结构网

 找回密码
 立即注册

QQ登录

只需一步,快速开始

膜结构车棚
膜结构车棚膜结构资质国产膜材 膜结构网中国膜结构协会
查看: 857|回复: 21

arx 点集计算

[复制链接]
  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

    发表于 2020-12-30 18:19 | 显示全部楼层 |阅读模式
    1. int pointSortByX(const void *p1, const void *p2){
    2.         return ((AcGePoint3d*)p1)->x - ((AcGePoint3d*)p2)->x;}

    3. int pointSortByY(const void *p1, const void *p2){
    4.         return ((AcGePoint3d*)p1)->y - ((AcGePoint3d*)p2)->y;}

    5. void sortPoints(AcGePoint3dArray allIntersPts)
    6. {
    7.         /*the first two elements are the start and endpoints. Use them to determine
    8.         if the array should be sorted by x or y values. */
    9.        
    10.         double distInX = abs(allIntersPts.at(0).x - allIntersPts.at(1).x);
    11.         double distInY = abs(allIntersPts.at(0).y - allIntersPts.at(1).y);
    12.        
    13.         //if the line is more vertical, then sort by the y-coords
    14.         if(distInY > distInX)
    15.                 std::qsort(allIntersPts.asArrayPtr(), allIntersPts.length(), sizeof(AcGePoint3d), pointSortByY);
    16.         //else the line is more horizontal, then sort by the x-coords
    17.         else
    18.                 std::qsort(allIntersPts.asArrayPtr(), allIntersPts.length(), sizeof(AcGePoint3d), pointSortByX);
    19.         return;
    20. }


    21. AcGePoint3dArray allIntersPts;
    22. /*code to append some points*/
    23. acutPrintf(_T("\nBefore sort:"));
    24.                 for(int i = 0; i < allIntersPts.logicalLength(); i++)
    25.                         acutPrintf(_T("\n%d: %.3f, %.3f"), i, allIntersPts.at(i).x, allIntersPts.at(i).y);

    26. sortPoints(allIntersPts);

    27. acutPrintf(_T("\nAfter sort:"));
    28.                 for(int i = 0; i < allIntersPts.logicalLength(); i++)
    29.                         acutPrintf(_T("\n%d: %.3f, %.3f"), i, allIntersPts.at(i).x, allIntersPts.at(i).y);
    复制代码
    回复


    http://www.mjgw.org/ 专业从事膜结构设计、制作加工、施工安装的膜结构工程服务,能够为客户提供专业的膜结构整体解决方案。做中国最好的膜结构综合服务平台。欢迎大家联系电话:198-7840-1958,QQ:463017170.
    相关关键词:膜结构车棚,膜结构车棚覆盖,膜结构车棚公司,膜结构车棚多少钱,膜结构车棚厂家,膜结构车棚价格,社区膜结构车棚,膜结构车棚膜布厂家 ,膜结构车棚哪家好,膜结构车棚多少钱一米,膜结构车棚报价,膜结构车棚哪里有,膜结构车棚定制,膜结构车棚安装,膜结构车棚设计,膜结构车棚电话,膜结构车棚加工,膜结构车棚膜布价格,膜结构车棚批发,膜结构车棚制造商,膜结构车棚生产厂家,膜结构车棚设计,膜结构车棚施工,膜结构车棚多少钱一平米,膜结构车棚订制,张拉膜车棚,张拉膜车棚覆盖,张拉膜车棚公司,张拉膜车棚多少钱,张拉膜车棚厂家,张拉膜车棚价格,社区张拉膜车棚,张拉膜车棚膜布厂家 ,张拉膜车棚哪家好,张拉膜车棚多少钱一米,张拉膜车棚报价,张拉膜车棚哪里有,张拉膜车棚定制,张拉膜车棚安装,张拉膜车棚设计,张拉膜车棚电话,张拉膜车棚加工,张拉膜车棚膜布价格,张拉膜车棚批发,张拉膜车棚制造商,张拉膜车棚生产厂家,张拉膜车棚设计,张拉膜车棚施工,张拉膜车棚多少钱一平米,张拉膜车棚订制,常用膜材品牌:德国杜肯、法国法拉利、德国海德斯、德国米乐、日本平岗、韩国秀博、比利时希运、美国赫虏伯、中国科宝、上海慧遥。

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2020-12-30 18:20 | 显示全部楼层
    1. std::sort(&data.first(), &data.last(), SortFileData);

    2. bool SorFileData(const PlotItem& item1, const PlotItem& item2)
    3. {
    4. double x1 = min(item1._start.x, item1._end.x);
    5. double x2 = min(item2._start.x, item2._end.x);

    6. return x1 < x2;
    7. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2020-12-30 18:20 | 显示全部楼层
    1. struct first_numeric_value
    2. {
    3. int ToInt(const std::string& str)
    4. {
    5. int x;
    6. std::stringstream ss(str);
    7. ss >> x;
    8. return x;
    9. }
    10. int getFirstNumericValue(const std::string& StrArray)
    11. {
    12. std::string::size_type pos1 = StrArray.find(",");
    13. std::string::size_type pos2 = StrArray.find(",",pos1+1);
    14. return ToInt(StrArray.substr(pos1+1,pos2-pos1-1));
    15. }

    16. bool operator()(const std::string& lhs, const std::string& rhs)
    17. {
    18. return getFirstNumericValue(lhs) < getFirstNumericValue(rhs);
    19. }
    20. };

    21. static void doSomeTest(void)
    22. {
    23. std::vector< std::string > myLst;
    24. std::vector< std::string >::iterator it;
    25. myLst.push_back( "\"T\",169.000,11.438,169.000,14.438,\"BO\"" );
    26. myLst.push_back( "\"P\",240.000,42.000,217.000,42.000,\"\"" );
    27. myLst.push_back( "\"T\",217.000,11.438,217.000,14.438,\"BO\"" );
    28. myLst.push_back( "\"P\",240.000,6.000,240.000,42.000,\"\"" );
    29. myLst.push_back( "\"T\",25.000,11.438,25.000,14.438,\"BO\"" );
    30. std::sort( myLst.begin(), myLst.end(), first_numeric_value() );
    31. for ( it = myLst.begin(); it != myLst.end(); ++it )
    32. {
    33. std::string s;
    34. s = ( *it );
    35. CString str( s.c_str() );
    36. acutPrintf( _T("\n%s"), str );
    37. }
    38. myLst.clear();
    39. }
    复制代码
    Then it will return an ascending sort of the vector:

    "T",25.000,11.438,25.000,14.438,"BO"
    "T",169.000,11.438,169.000,14.438,"BO"
    "T",217.000,11.438,217.000,14.438,"BO"
    "P",240.000,42.000,217.000,42.000,""
    "P",240.000,6.000,240.000,42.000,""[/code]
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2020-12-30 18:21 | 显示全部楼层
    1. PointSorter(const AcGePoint3d& anchor, AcGePoint3dArray& points)
    2. {
    3.         mAnchor = anchor;
    4.         AcGePoint3d* arr = points.asArrayPtr();
    5.         //auto size = points.length() / sizeof(AcGePoint3d); //WRONG!
    6.         auto size = points.length();

    7.         std::sort(arr, arr + size,
    8.                 [this](const AcGePoint3d& l, const AcGePoint3d& r)
    9.                 -> bool { return (doCompare(l, r) < 0); } // true if l<r,  false if l>=r
    10.         );
    11. }

    12. PointSorter(const AcGeVector3d& axis, const AcGePoint3d& anchor, AcGePoint3dArray& points)
    13. {
    14.         mAxis = axis;
    15.         mAnchor = anchor;
    16.         AcGePoint3d* arr = points.asArrayPtr();
    17.         // auto size = points.length() / sizeof(AcGePoint3d); //WRONG!
    18.         auto size = points.length();

    19.         std::sort(arr, arr + size,
    20.                 [this](const AcGePoint3d& l, const AcGePoint3d& r)
    21.                 -> bool { return (doCompareAlongAxis(l, r) < 0); } // true if l<r,  false if l>=r
    22.         );
    23. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 09:57 | 显示全部楼层
    1. //返回 点与点的平面距离
    2. double TwoDistancePointAndPoint(const xjPoint &p1, const xjPoint &p2)
    3. {
    4.         double x1 = p1.x, y1 = p1.y;
    5.         double x2 = p2.x, y2 = p2.y;

    6.         double dis = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));

    7.         return dis;
    8. }

    9. //返回 点(p0)  到线(p1,p2)的距离
    10. double TwoDistancePointAndLine(const xjPoint &p0, const xjPoint &p1,const xjPoint &p2)
    11. {
    12.         double dis12 = TwoDistancePointAndPoint(p1, p2);//线段长度
    13.         double dis01 = TwoDistancePointAndPoint(p0, p1);//p1与p0的距离
    14.         double dis02 = TwoDistancePointAndPoint(p0, p2);//p2与p0的距离
    15.         double HalfC = (dis12 + dis01 + dis02) / 2;// 半周长
    16.         double s = sqrt(HalfC * (HalfC - dis12) * (HalfC - dis01) * (HalfC - dis02));//海伦公式求面积
    17.         double xj2DisPL = 2 * s / dis12;// 返回点到线的距离(利用三角形面积公式求高)

    18.         return xj2DisPL;
    19. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 09:58 | 显示全部楼层
    1. /*返回 最小外接矩形点*/
    2. /* xjListCH 凸包点 */
    3. void xjMinimumBoundingRectangle(QList<xjPoint> &xjListMBRpoint, double &length, double &width, QList<xjPoint> xjListCH)
    4. {
    5.         //外接矩形
    6.         QList<xjPoint> xjListCH2;
    7.         int xjPcount = xjListCH.size();
    8.         for (int i = 0; i < xjPcount; i++)
    9.         {
    10.                 xjListCH2.append(xjListCH.at(i));
    11.         }
    12.         qSort(xjListCH2.begin(), xjListCH2.end(), [](const xjPoint &a, const xjPoint &b) {return a.x < b.x; });
    13.         double minX = xjListCH2.at(0).x;
    14.         double maxX = xjListCH2.at(xjPcount - 1).x;
    15.         qSort(xjListCH2.begin(), xjListCH2.end(), [](const xjPoint &a, const xjPoint &b) {return a.y < b.y; });
    16.         double minY = xjListCH2.at(0).y;
    17.         double maxY = xjListCH2.at(xjPcount - 1).y;

    18.         //依次判断
    19.         double minArea = 99999999;
    20.         xjListCH.push_back(xjListCH.at(0));
    21.         for (int a = 0; a < xjListCH.size() - 1; a++)
    22.         {
    23.                 xjPoint p0 = xjListCH.at(a);
    24.                 xjPoint p1 = xjListCH.at(a + 1);

    25.                 if ((p0.y == p1.y)|| (p0.x == p1.x)) // 水平或垂直
    26.                 {
    27.                         double side1 = maxY - minY;
    28.                         double side0 = maxX - minX;
    29.                         double xjArea = side0 * side1;

    30.                         if (xjArea <= minArea)
    31.                         {
    32.                                 length = max(side0, side1);
    33.                                 width = min(side0, side1);
    34.                                 minArea = xjArea;

    35.                                 //外接矩形四个点
    36.                                 xjPoint pLB;
    37.                                 pLB.x = minX;
    38.                                 pLB.y = minY;
    39.                                 pLB.z = 0;
    40.                                 xjPoint pRB;
    41.                                 pRB.x = maxX;
    42.                                 pRB.y = minY;
    43.                                 pRB.z = 0;
    44.                                 xjPoint pRT;
    45.                                 pRT.x = maxX;
    46.                                 pRT.y = maxY;
    47.                                 pRT.z = 0;
    48.                                 xjPoint pLT;
    49.                                 pLT.x = minX;
    50.                                 pLT.y = maxY;
    51.                                 pLT.z = 0;

    52.                                 xjListMBRpoint.clear();
    53.                                 xjListMBRpoint.append(pLB);
    54.                                 xjListMBRpoint.append(pRB);
    55.                                 xjListMBRpoint.append(pRT);
    56.                                 xjListMBRpoint.append(pLT);
    57.                         }
    58.                 }
    59.                 else //不水平 不垂直
    60.                 {
    61.                         double k1 = (p1.y - p0.y) / (p1.x - p0.x);
    62.                         double b1 = p0.y - k1 * p0.x;
    63.                         double side0 = -3;
    64.                         xjPoint Pside0;
    65.                         for (int j = 0; j < xjListCH.size(); j++)
    66.                         {
    67.                                 if ((j == a) || (j == (a + 1)))
    68.                                         continue;

    69.                                 xjPoint p = xjListCH.at(j);
    70.                                 double dis = abs(TwoDistancePointAndLine(p, p0, p1));
    71.                                 if (dis >= side0)
    72.                                 {
    73.                                         side0 = dis;
    74.                                         Pside0.x = p.x;
    75.                                         Pside0.y = p.y;
    76.                                         Pside0.z = p.z;
    77.                                 }
    78.                         }
    79.                         double b11 = Pside0.y - k1 * Pside0.x;

    80.                         //垂直方向
    81.                         double k2 = -1.0 / k1;
    82.                         double bb = p0.y - k2 * p0.x;
    83.                         double side1_positive = -3;
    84.                         xjPoint Pside1_positive;
    85.                         double side1_negative = 9999999;
    86.                         xjPoint Pside1_negative;
    87.                         for (int j = 0; j < xjListCH.size(); j++)
    88.                         {
    89.                                 xjPoint p = xjListCH.at(j);
    90.                                 double dis = (k2*p.x - p.y + bb) / (sqrt(k2*k2 + 1));
    91.                                 if ((dis>=0)&&(dis >= side1_positive))
    92.                                 {
    93.                                         side1_positive = dis;
    94.                                         Pside1_positive.x = p.x;
    95.                                         Pside1_positive.y = p.y;
    96.                                         Pside1_positive.z = p.z;
    97.                                 }
    98.                                 if ((dis<0)&&(dis <= side1_negative))
    99.                                 {
    100.                                         side1_negative = dis;
    101.                                         Pside1_negative.x = p.x;
    102.                                         Pside1_negative.y = p.y;
    103.                                         Pside1_negative.z = p.z;
    104.                                 }
    105.                         }

    106.                         double b2 = Pside1_positive.y - k2 * Pside1_positive.x;
    107.                         double b22 = Pside1_negative.y - k2 * Pside1_negative.x;

    108.                         //面积和周长
    109.                         double side1 = abs(side1_positive)+abs(side1_negative);
    110.                         double xjArea = side0 * side1;
    111.                        
    112.                         if (xjArea <= minArea)
    113.                         {
    114.                                 length = max(side0, side1);
    115.                                 width = min(side0, side1);
    116.                                 minArea = xjArea;

    117.                                 //外接矩形四个点
    118.                                 xjPoint br0;
    119.                                 br0.x = (b1 - b22) / (k2 - k1);
    120.                                 br0.y = k1 * br0.x + b1;
    121.                                 xjPoint br1;
    122.                                 br1.x = (b11-b22) / (k2-k1);
    123.                                 br1.y =  k1* br1.x + b11;
    124.                                 xjPoint br2;
    125.                                 br2.x = (b2-b11) / (k1-k2);
    126.                                 br2.y =  k1* br2.x + b11;
    127.                                 xjPoint br3;
    128.                                 br3.x = (b2-b1) / (k1-k2);
    129.                                 br3.y = k1 * br3.x + b1;

    130.                                 xjListMBRpoint.clear();
    131.                                 xjListMBRpoint.append(br0);
    132.                                 xjListMBRpoint.append(br1);
    133.                                 xjListMBRpoint.append(br2);
    134.                                 xjListMBRpoint.append(br3);
    135.                         }
    136.                 }
    137.         }

    138.         //MBR提示信息
    139.         QString MBRinfo = "chMBR: length = " + QString::number(length, 'f', 4) + ", ";
    140.         MBRinfo += "width = " + QString::number(width, 'f', 4) + ", ";
    141.         MBRinfo += "minimum area = " + QString::number(minArea) + ", ";
    142.         MBRinfo += "circumference = " + QString::number((length + width) * 2);
    143. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 09:59 | 显示全部楼层
    1. //自定义结构体
    2. struct xjPoint
    3. {
    4.         double x;
    5.         double y;
    6.         double z;
    7. };


    8. /*返回 两点的平面距离*/
    9. double xjGet2DistancePP(xjPoint p1, xjPoint p2)
    10. {
    11.         return (sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)));
    12. }


    13. /*返回 两个向量的叉积(p0是公共点)*/
    14. double xjGetCrossProduct(xjPoint p1, xjPoint p2, xjPoint p0)
    15. {
    16.         return ((p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y));
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 10:00 | 显示全部楼层
    1. /*方式一、返回 凸包点集,通过 Graham Scan 方法*/
    2. QList<xjPoint> xjGetConvexHullByGrahamScan(QList<xjPoint> listPoint)
    3. {
    4.         QList<xjPoint> xjListResult;

    5.         int i = 0, j = 0, k = 0;
    6.         xjPoint tmP = listPoint[0];

    7.         //排序 找到最下且偏左的点
    8.         for (i = 1; i < listPoint.size(); i++)
    9.         {
    10.                 if ((listPoint[i].y < tmP.y) || ((listPoint[i].y == tmP.y) && (listPoint[i].x < tmP.x)))
    11.                 {
    12.                         tmP = listPoint[i];
    13.                         k = i;
    14.                 }
    15.         }
    16.         listPoint[k] = listPoint[0];
    17.         listPoint[0] = tmP;

    18.         //排序 按极角从小到大,距离从近到远
    19.         for (i = 1; i < listPoint.size(); i++)
    20.         {
    21.                 k = i;
    22.                 for (j = i + 1; j < listPoint.size(); j++)
    23.                 {
    24.                         double cross = xjGetCrossProduct(listPoint[j], listPoint[k], listPoint[0]);
    25.                         double disance_j = xjGet2DistancePP(listPoint[0], listPoint[j]);
    26.                         double disance_k = xjGet2DistancePP(listPoint[0], listPoint[k]);
    27.                         if ((cross > 0) || ((cross == 0) && (disance_j < disance_k)))
    28.                         {
    29.                                 k = j;//k保存极角最小的那个点,或者相同距离原点最近
    30.                         }
    31.                 }
    32.                 tmP = listPoint[i];
    33.                 listPoint[i] = listPoint[k];
    34.                 listPoint[k] = tmP;
    35.         }

    36.         //添加第一、二个凸包点
    37.         xjListResult.append(listPoint[0]);
    38.         xjListResult.append(listPoint[1]);

    39.         遍历
    40.         for (i = 2; i < listPoint.size(); ++i)
    41.         {
    42.                 for (j = i + 1; j < listPoint.size(); ++j)
    43.                 {
    44.                         if (xjGetCrossProduct(listPoint[i], listPoint[j], xjListResult[xjListResult.size() - 1]) < 0)
    45.                         {
    46.                                 i = j;
    47.                         }
    48.                 }
    49.                 xjListResult.append(listPoint[i]);
    50.         }

    51.         return xjListResult;
    52. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 10:00 | 显示全部楼层
    1. /*方式二:返回 凸包点集,通过 Graham Scan 方法*/
    2. QList<xjPoint> xjGetConvexHull(QList<xjPoint> list)
    3. {
    4.          QList<xjPoint> newList;

    5.          /*最下偏左的点 肯定是凸包点*/
    6.          qSort(list.begin(), list.end(), [](const xjPoint &a, const xjPoint &b)
    7.          {
    8.                  if (a.y == b.y)
    9.                          return a.x < b.x;
    10.                  else
    11.                          return a.y < b.y;
    12.          });
    13.          xjPoint p0 = list.at(0);

    14.          //临时点集 坐标偏移
    15.          QList<xjPoint> tempList;
    16.          for (int i = 0; i < list.size(); i++)
    17.          {
    18.                  xjPoint p;
    19.                  p.x = list.at(i).x - p0.x;
    20.                  p.y = list.at(i).y - p0.y;
    21.                  p.z = list.at(i).z - p0.z;
    22.                  tempList.append(p);
    23.          }
    24.          newList.append(tempList.at(0));

    25.          //按极角从小到大,距离偏短进行排序
    26.          qSort(tempList.begin(), tempList.end(), [](const xjPoint &a, const xjPoint &b)
    27.          {
    28.                  if (atan2(a.y, a.x) != atan2(b.y, b.x))
    29.                          return (atan2(a.y, a.x)) < (atan2(b.y, b.x));
    30.                  return a.x < b.x;
    31.          });
    32.          xjPoint p1 = tempList.at(1);
    33.          newList.append(p1);

    34.          //获取凸包点
    35.          for (int j = 2; j < tempList.size(); j++)
    36.          {
    37.                  for (int k = j + 1; k < tempList.size(); k++)
    38.                  {
    39.                          double xjCross = xjGetCrossProduct(tempList.at(j), tempList.at(k), newList.at(newList.size() - 1));
    40.                          if (xjCross < 0)
    41.                          {
    42.                                  j = k;
    43.                          }
    44.                  }
    45.                  newList.append(tempList.at(j));
    46.          }


    47.          //回复坐标
    48.          for (int i = 0; i < newList.size(); i++)
    49.          {
    50.                  xjPoint p = newList.at(i);
    51.                  p.x += p0.x;
    52.                  p.y += p0.y;
    53.                  p.z += p0.z;
    54.                  newList.replace(i, p);
    55.          }

    56.          return newList;
    57. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 15:34 | 显示全部楼层
    1. // 判断两条多段线是否重叠//
    2. // 返回值:-1——(pl1<pl2)
    3. //                  0——(pl1=pl2)
    4. //                    1——(pl1>pl2)
    5. int PlineOverlap(AcDbPolyline* pl1, AcDbPolyline* pl2)
    6. {
    7.         int rtn1(1), rtn2(1);
    8.         Acad::ErrorStatus es1, es2;
    9.         AcGePoint3d pt1, pt2;
    10.         double param1(0), param2(0);      
    11.         int i,j;      
    12.         // 先判断其中一条线的所有结点是否在另一条线上
    13.         for (i=0,j=1; j<pl1->numVerts(); i++,j++)
    14.         {
    15.                 pl1->getPointAt(i, pt1);
    16.                 pl1->getPointAt(j, pt2);
    17.                 if (pl2->getParamAtPoint(pt1, param1)!=Acad::eOk
    18.                         || pl2->getParamAtPoint(pt2, param2)!=Acad::eOk)
    19.                 {
    20.                         es1 = Acad::ePointNotOnEntity;
    21.                         break;
    22.                 }
    23.                 else
    24.                 {        //再判断两条线的子段是否有交点                                                                              
    25.                         AcDbPolyline* pl2Sub = GetSplitPline(pl2, param1, param2);
    26.                         if (HasInter(pt1, pt2, pl2Sub))
    27.                         {
    28.                                 break;
    29.                         }                                               
    30.                         rtn1 ++;                       
    31.                 }
    32.         }
    33.         for (i=0,j=1; j<pl2->numVerts(); i++,j++)
    34.         {
    35.                 pl2->getPointAt(i, pt1);
    36.                 pl2->getPointAt(j, pt2);
    37.                 if (pl1->getParamAtPoint(pt1, param1)!=Acad::eOk
    38.                         || pl1->getParamAtPoint(pt2, param2)!=Acad::eOk)
    39.                 {
    40.                         es2 = Acad::ePointNotOnEntity;
    41.                         break;
    42.                 }
    43.                 else
    44.                 {                                                      
    45.                         AcDbPolyline* pl1Sub = GetSplitPline(pl1, param1, param2);
    46.                         if (HasInter(pt1, pt2, pl1Sub))
    47.                         {
    48.                                 break;
    49.                         }                       
    50.                         rtn2 ++;                                               
    51.                 }
    52.         }
    53.                        
    54.         if (rtn1==pl1->numVerts()
    55.                 && rtn2==pl2->numVerts())
    56.         {
    57.                 return 0;
    58.         }
    59.         else if (rtn1==pl1->numVerts())
    60.         {
    61.                 return -1;
    62.         }
    63.         else if (rtn2==pl2->numVerts())
    64.         {
    65.                 return 1;
    66.         }
    67.         else
    68.         {
    69.                 return (-99);
    70.         }                       
    71. }
    72. //判断线段与多段线是否有交点
    73. BOOL HasInter(AcGePoint3d ptFrom, AcGePoint3d ptTo, AcDbPolyline* pl, BOOL bExtend=FALSE)
    74. {
    75.         int i,j;
    76.         AcGePoint3d pt1, pt2;
    77.         for (i=0,j=1; j<pl->numVerts(); i++,j++)
    78.         {
    79.                 pl->getPointAt(i, pt1);
    80.                 pl->getPointAt(j, pt2);
    81.                 ads_point inter;
    82.                 int teston;
    83.                 if (bExtend) teston=0;
    84.                 else teston=1;
    85.                 if (acdbInters(asDblArray(ptFrom), asDblArray(ptTo), asDblArray(pt1), asDblArray(pt2), teston, inter)==RTNORM)
    86.                 {
    87.                         return TRUE;
    88.                 }
    89.         }
    90.         return FALSE;
    91. }

    92. //取中间段
    93. AcDbPolyline* GetSplitPline(AcDbPolyline* pl, double param1, double param2)
    94. {
    95.         AcGeDoubleArray params;
    96.         if (param1<param2)
    97.         {
    98.                 params.append(param1);
    99.                 params.append(param2);
    100.         }
    101.         else
    102.         {
    103.                 params.append(param2);
    104.                 params.append(param1);
    105.         }
    106.         AcDbVoidPtrArray segs;
    107.         if (pl->getSplitCurves(params, segs)==Acad::eOk)
    108.         {
    109.                 if (segs.length()==1)
    110.                 {
    111.                         return static_cast<AcDbPolyline*>(segs.at(0));
    112.                 }
    113.                 else if (segs.length()==2)
    114.                 {
    115.                         if (params.at(0)==0)
    116.                         {
    117.                                 return static_cast<AcDbPolyline*>(segs.at(0));
    118.                         }
    119.                         if (params.at(1)==(pl->numVerts()-1))
    120.                         {
    121.                                 return static_cast<AcDbPolyline*>(segs.at(1));
    122.                         }
    123.                 }
    124.                 else if (segs.length()==3)
    125.                 {
    126.                         return static_cast<AcDbPolyline*>(segs.at(1));
    127.                 }
    128.         }
    129.         return pl;
    130. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-10 15:34 | 显示全部楼层
    1. 用的时候需要加入头文件#include <math.h>
    2. 或者#include <math>
    3. using namespace std;
    4. 常用函数:

    5. abs绝对值函数
    6. acos反余弦函数
    7. asin反正弦函数
    8. atan反正切函数
    9. exp e的x次方
    10. cos余弦函数
    11. sin正弦函数
    12. tan正切函数
    13. ceil求不小于x的最小整数
    14. cosh求x的双曲余弦值
    15. fabs求浮点数x的绝对值
    16. fmod计算x/y的余数
    17. frexp把浮点数x分解成尾数和指数
    18. hypot对于给定的直角三角形的两个直角边,求其斜边的长度。
    19. ldexp装载浮点数
    20. log,e为底对数
    21. log10,10为底对数
    22. modf,把数分为指数和尾数
    23. pow,计算x的y次幂
    24. sinh,求x的双曲正弦值
    25. sqrt,开方
    26. tanh求x的双曲正切值
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-14 16:24 | 显示全部楼层
    1. //点表比较函数,配合std::sort使用
    2. //点按xyz排序,小写字母表示从小到大排序,大写则从大到小排序,在前的字母先排序
    3. int cmp_xyz (AcGePoint3d p1,AcGePoint3d p2)
    4. {
    5.         if (p1.x != p2.x)
    6.                 return p1.x < p2.x?1:0;
    7.         else if (p1.y != p2.y)
    8.                 return p1.y < p2.y?1:0;
    9.         else
    10.                 return p1.z < p2.z?1:0;
    11. }
    12. int cmp_xYz (AcGePoint3d p1,AcGePoint3d p2)
    13. {
    14.         if (p1.x != p2.x)
    15.                 return p1.x < p2.x?1:0;
    16.         else if (p1.y != p2.y)
    17.                 return p1.y > p2.y?1:0;
    18.         else
    19.                 return p1.z < p2.z?1:0;
    20. }
    21. int cmp_xYZ (AcGePoint3d p1,AcGePoint3d p2)
    22. {
    23.         if (p1.x != p2.x)
    24.                 return p1.x < p2.x?1:0;
    25.         else if (p1.y != p2.y)
    26.                 return p1.y > p2.y?1:0;
    27.         else
    28.                 return p1.z > p2.z?1:0;
    29. }
    30. int cmp_xyZ (AcGePoint3d p1,AcGePoint3d p2)
    31. {
    32.         if (p1.x != p2.x)
    33.                 return p1.x < p2.x?1:0;
    34.         else if (p1.y != p2.y)
    35.                 return p1.y < p2.y?1:0;
    36.         else
    37.                 return p1.z > p2.z?1:0;
    38. }
    39. int cmp_XYZ (AcGePoint3d p1,AcGePoint3d p2)
    40. {
    41.         if (p1.x != p2.x)
    42.                 return p1.x > p2.x?1:0;
    43.         else if (p1.y != p2.y)
    44.                 return p1.y > p2.y?1:0;
    45.         else
    46.                 return p1.z > p2.z?1:0;
    47. }
    48. int cmp_XyZ (AcGePoint3d p1,AcGePoint3d p2)
    49. {
    50.         if (p1.x != p2.x)
    51.                 return p1.x > p2.x?1:0;
    52.         else if (p1.y != p2.y)
    53.                 return p1.y < p2.y?1:0;
    54.         else
    55.                 return p1.z > p2.z?1:0;
    56. }

    57. int cmp_Xyz (AcGePoint3d p1,AcGePoint3d p2)
    58. {
    59.         if (p1.x != p2.x)
    60.                 return p1.x > p2.x?1:0;
    61.         else if (p1.y != p2.y)
    62.                 return p1.y < p2.y?1:0;
    63.         else
    64.                 return p1.z < p2.z?1:0;
    65. }

    66. int cmp_XYz (AcGePoint3d p1,AcGePoint3d p2)
    67. {
    68.         if (p1.x != p2.x)
    69.                 return p1.x > p2.x?1:0;
    70.         else if (p1.y != p2.y)
    71.                 return p1.y > p2.y?1:0;
    72.         else
    73.                 return p1.z < p2.z?1:0;
    74. }

    75. int cmp_xy (AcGePoint3d p1,AcGePoint3d p2)
    76. {
    77.         if (p1.x != p2.x)
    78.                 return p1.x < p2.x?1:0;
    79.         else
    80.                 return p1.y < p2.y?1:0;        
    81. }
    82. int cmp_xY (AcGePoint3d p1,AcGePoint3d p2)
    83. {
    84.         if (p1.x != p2.x)
    85.                 return p1.x < p2.x?1:0;
    86.         else
    87.                 return p1.y > p2.y?1:0;        
    88. }
    89. int cmp_Xy (AcGePoint3d p1,AcGePoint3d p2)
    90. {
    91.         if (p1.x != p2.x)
    92.                 return p1.x > p2.x?1:0;
    93.         else
    94.                 return p1.y < p2.y?1:0;        
    95. }
    96. int cmp_XY (AcGePoint3d p1,AcGePoint3d p2)
    97. {
    98.         if (p1.x != p2.x)
    99.                 return p1.x > p2.x?1:0;
    100.         else
    101.                 return p1.y > p2.y?1:0;        
    102. }

    103. int cmp_x(AcGePoint3d p1,AcGePoint3d p2)
    104. {
    105.         return p1.x < p2.x?1:0;
    106. }
    107. int cmp_X(AcGePoint3d p1,AcGePoint3d p2)
    108. {
    109.         return p1.x > p2.x?1:0;
    110. }
    111. int cmp_y(AcGePoint3d p1,AcGePoint3d p2)
    112. {
    113.         return p1.y < p2.y?1:0;
    114. }
    115. int cmp_Y(AcGePoint3d p1,AcGePoint3d p2)
    116. {
    117.         return p1.y > p2.y?1:0;
    118. }
    119. int cmp_z(AcGePoint3d p1,AcGePoint3d p2)
    120. {
    121.         return p1.z < p2.z?1:0;
    122. }
    123. int cmp_Z(AcGePoint3d p1,AcGePoint3d p2)
    124. {
    125.         return p1.z > p2.z?1:0;
    126. }
    127. int cmp_YX (AcGePoint3d p1,AcGePoint3d p2)
    128. {
    129.         if (p1.y != p2.y)
    130.                 return p1.y > p2.y?1:0;
    131.         else
    132.                 return p1.x > p2.x?1:0;        
    133. }
    134. int cmp_Yx (AcGePoint3d p1,AcGePoint3d p2)
    135. {
    136.         if (p1.y != p2.y)
    137.                 return p1.y > p2.y?1:0;
    138.         else
    139.                 return p1.x < p2.x?1:0;        
    140. }
    141. int cmp_yx (AcGePoint3d p1,AcGePoint3d p2)
    142. {
    143.         if (p1.y != p2.y)
    144.                 return p1.y < p2.y?1:0;
    145.         else
    146.                 return p1.x < p2.x?1:0;        
    147. }
    148. int cmp_yX (AcGePoint3d p1,AcGePoint3d p2)
    149. {
    150.         if (p1.y != p2.y)
    151.                 return p1.y < p2.y?1:0;
    152.         else
    153.                 return p1.x > p2.x?1:0;        
    154. }
    155. int cmp_yxz (AcGePoint3d p1,AcGePoint3d p2)
    156. {
    157.         if (p1.y != p2.y)
    158.                 return p1.y < p2.y?1:0;
    159.         else if (p1.x != p2.x)
    160.                 return p1.x < p2.x?1:0;
    161.         else
    162.                 return p1.z < p2.z?1:0;
    163. }
    164. int cmp_yXz (AcGePoint3d p1,AcGePoint3d p2)
    165. {
    166.         if (p1.y != p2.y)
    167.                 return p1.y < p2.y?1:0;
    168.         else if (p1.x != p2.x)
    169.                 return p1.x > p2.x?1:0;
    170.         else
    171.                 return p1.z < p2.z?1:0;
    172. }
    173. int cmp_yXZ (AcGePoint3d p1,AcGePoint3d p2)
    174. {
    175.         if (p1.y != p2.y)
    176.                 return p1.y < p2.y?1:0;
    177.         else if (p1.x != p2.x)
    178.                 return p1.x > p2.x?1:0;
    179.         else
    180.                 return p1.z > p2.z?1:0;
    181. }
    182. int cmp_yxZ (AcGePoint3d p1,AcGePoint3d p2)
    183. {
    184.         if (p1.y != p2.y)
    185.                 return p1.y < p2.y?1:0;
    186.         else if (p1.x != p2.x)
    187.                 return p1.x < p2.x?1:0;
    188.         else
    189.                 return p1.z > p2.z?1:0;
    190. }
    191. int cmp_YXZ (AcGePoint3d p1,AcGePoint3d p2)
    192. {
    193.         if (p1.y != p2.y)
    194.                 return p1.y > p2.y?1:0;
    195.         else if (p1.x != p2.x)
    196.                 return p1.x > p2.x?1:0;
    197.         else
    198.                 return p1.z > p2.z?1:0;
    199. }
    200. int cmp_YxZ (AcGePoint3d p1,AcGePoint3d p2)
    201. {
    202.         if (p1.y != p2.y)
    203.                 return p1.y > p2.y?1:0;
    204.         else if (p1.x != p2.x)
    205.                 return p1.x < p2.x?1:0;
    206.         else
    207.                 return p1.z > p2.z?1:0;
    208. }

    209. int cmp_Yxz (AcGePoint3d p1,AcGePoint3d p2)
    210. {
    211.         if (p1.y != p2.y)
    212.                 return p1.y > p2.y?1:0;
    213.         else if (p1.x != p2.x)
    214.                 return p1.x < p2.x?1:0;
    215.         else
    216.                 return p1.z < p2.z?1:0;
    217. }

    218. int cmp_YXz (AcGePoint3d p1,AcGePoint3d p2)
    219. {
    220.         if (p1.y != p2.y)
    221.                 return p1.y > p2.y?1:0;
    222.         else if (p1.x != p2.x)
    223.                 return p1.x > p2.x?1:0;
    224.         else
    225.                 return p1.z < p2.z?1:0;
    226. }


    227. int cmp_2d_xy (AcGePoint2d p1,AcGePoint2d p2)
    228. {
    229.         if (p1.x != p2.x)
    230.                 return p1.x < p2.x?1:0;
    231.         else
    232.                 return p1.y < p2.y?1:0;        
    233. }
    234. int cmp_2d_xY (AcGePoint2d p1,AcGePoint2d p2)
    235. {
    236.         if (p1.x != p2.x)
    237.                 return p1.x < p2.x?1:0;
    238.         else
    239.                 return p1.y > p2.y?1:0;        
    240. }
    241. int cmp_2d_Xy (AcGePoint2d p1,AcGePoint2d p2)
    242. {
    243.         if (p1.x != p2.x)
    244.                 return p1.x > p2.x?1:0;
    245.         else
    246.                 return p1.y < p2.y?1:0;        
    247. }
    248. int cmp_2d_XY (AcGePoint2d p1,AcGePoint2d p2)
    249. {
    250.         if (p1.x != p2.x)
    251.                 return p1.x > p2.x?1:0;
    252.         else
    253.                 return p1.y > p2.y?1:0;        
    254. }

    255. int cmp_2d_x (AcGePoint2d p1,AcGePoint2d p2)
    256. {
    257.         return p1.x < p2.x?1:0;
    258. }
    259. int cmp_2d_X(AcGePoint2d p1,AcGePoint2d p2)
    260. {
    261.         return p1.x > p2.x?1:0;
    262. }
    263. int cmp_2d_y (AcGePoint2d p1,AcGePoint2d p2)
    264. {
    265.         return p1.y < p2.y?1:0;
    266. }
    267. int cmp_2d_Y(AcGePoint2d p1,AcGePoint2d p2)
    268. {
    269.         return p1.y > p2.y?1:0;
    270. }

    271. int cmp_2d_YX(AcGePoint2d p1,AcGePoint2d p2)
    272. {
    273.         if (p1.y != p2.y)
    274.                 return p1.y > p2.y?1:0;
    275.         else
    276.                 return p1.x > p2.x?1:0;        
    277. }
    278. int cmp_2d_Yx (AcGePoint2d p1,AcGePoint2d p2)
    279. {
    280.         if (p1.y != p2.y)
    281.                 return p1.y > p2.y?1:0;
    282.         else
    283.                 return p1.x < p2.x?1:0;        
    284. }
    285. int cmp_2d_yx (AcGePoint2d p1,AcGePoint2d p2)
    286. {
    287.         if (p1.y != p2.y)
    288.                 return p1.y < p2.y?1:0;
    289.         else
    290.                 return p1.x < p2.x?1:0;        
    291. }
    292. int cmp_2d_yX (AcGePoint2d p1,AcGePoint2d p2)
    293. {
    294.         if (p1.y != p2.y)
    295.                 return p1.y < p2.y?1:0;
    296.         else
    297.                 return p1.x > p2.x?1:0;        
    298. }

    299. int cmp_up (int &a,int &b)
    300. {
    301.         return a<b?1:0;
    302. }
    303. int cmp_up (double &a,double &b)
    304. {
    305.         return a<b?1:0;
    306. }
    307. int cmp_down (const int &a, const int &b)
    308. {
    309.         return a>b?1:0;
    310. }
    311. int cmp_down (double &a,double &b)
    312. {
    313.         return a>b?1:0;
    314. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-14 16:25 | 显示全部楼层
    1. int cmp_2Point (AcGePoint3d p1,AcGePoint3d p2,CString &cmpstr)
    2. {
    3.   if (p1.x != p2.x)
    4.     return  acutWcMatch(cmpstr, _T("X*") )== RTNORM? p1.x > p2.x : p1.x < p2.x ;
    5.   else if (p1.y != p2.y)
    6.       return  acutWcMatch(cmpstr, _T("?Y*") )== RTNORM? p1.y > p2.y : p1.y < p2.y;
    7.   else
    8.     return  acutWcMatch(cmpstr, _T("*Z") )== RTNORM? p1.z > p2.z : p1.z < p2.z;
    9. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-14 16:26 | 显示全部楼层
    1. /*******************************************************************************
    2. * Filename:    multiway_qsort.hpp
    3. * Author:      Zhang Li
    4. * Descrption:
    5. *      Faster soring a set of strings with long common prefix (like URLs).
    6. *
    7. ******************************************************************************/
    8. #ifndef HEADER_MULTIWAY_QSORT_HPP
    9. #define HEADER_MULTIWAY_QSORT_HPP

    10. #include <iterator>

    11. template<class RandIt> static inline void multiway_qsort(RandIt beg, RandIt end, size_t depth = 0, size_t step_len = 6) {
    12.     if(beg + 1 >= end) {
    13.         return;
    14.     }

    15.     struct { /* implement bounded comparing */
    16.         inline int operator() (
    17.                 const typename std::iterator_traits<RandIt>::value_type& a,
    18.                 const typename std::iterator_traits<RandIt>::value_type& b, size_t depth, size_t step_len) {

    19.             for(size_t i = 0; i < step_len; i++) {
    20.                 if(a[depth + i] == b[depth + i] && a[depth + i] == 0) return 0;
    21.                 if(a[depth + i] <  b[depth + i]) return +1;
    22.                 if(a[depth + i] >  b[depth + i]) return -1;
    23.             }
    24.             return 0;
    25.         }
    26.     } bounded_cmp;

    27.     RandIt i = beg;
    28.     RandIt j = beg + std::distance(beg, end) / 2;
    29.     RandIt k = end - 1;

    30.     typename std::iterator_traits<RandIt>::value_type key = ( /* median of l,m,r */
    31.             bounded_cmp(*i, *j, depth, step_len) > 0 ?
    32.             (bounded_cmp(*i, *k, depth, step_len) > 0 ? (bounded_cmp(*j, *k, depth, step_len) > 0 ? *j : *k) : *i) :
    33.             (bounded_cmp(*i, *k, depth, step_len) < 0 ? (bounded_cmp(*j, *k, depth, step_len) < 0 ? *j : *k) : *i));

    34.     /* 3-way partition */
    35.     for(j = i; j <= k; ++j) {
    36.         switch(bounded_cmp(*j, key, depth, step_len)) {
    37.             case +1: std::iter_swap(i, j); ++i;      break;
    38.             case -1: std::iter_swap(k, j); --k; --j; break;
    39.         }
    40.     }
    41.     ++k;

    42.     if(beg + 1 < i) multiway_qsort(beg, i, depth, step_len); /* recursively sort [x > pivot] subset */
    43.     if(end + 1 > k) multiway_qsort(k, end, depth, step_len); /* recursively sort [x < pivot] subset */

    44.     /* recursively sort [x == pivot] subset with higher depth */
    45.     if(i < k && (*i)[depth] != 0) {
    46.         multiway_qsort(i, k, depth + step_len, step_len);
    47.     }
    48.     return;
    49. }
    50. #endif

    51. #include <boost/progress.hpp>
    52. #include <boost/utility/string_ref.hpp>
    53. #include <iostream>
    54. #include <vector>
    55. #include <cassert>

    56. int main() {
    57.     std::vector<std::string> items;
    58.     std::string s;
    59.     std::ios::sync_with_stdio(false);

    60.     /* init items from input */
    61.     while(std::getline(std::cin, s)) {
    62.         items.push_back(s);
    63.     }

    64.     std::vector<boost::string_ref> refs_1(items.begin(), items.end());
    65.     std::vector<boost::string_ref> refs_2(items.begin(), items.end());
    66.     std::vector<boost::string_ref> refs_3(items.begin(), items.end());

    67.     /* multiway_qsort */
    68.     std::cerr << "sorting with multiway_qsort()...\n";
    69.     {
    70.         boost::progress_timer _(std::cerr);
    71.         multiway_qsort(refs_1.begin(), refs_1.end());
    72.     }

    73.     /* std::sort */
    74.     std::cerr << "sorting with std::sort()...\n";
    75.     {
    76.         boost::progress_timer _(std::cerr);
    77.         std::sort(refs_2.begin(), refs_2.end());
    78.     }

    79.     /* std::stable_sort */
    80.     std::cerr << "sorting with std::stable_sort()...\n";
    81.     {
    82.         boost::progress_timer _(std::cerr);
    83.         std::stable_sort(refs_3.begin(), refs_3.end());
    84.     }

    85.     /* verify sorting result of multiway_sort */
    86.     assert(refs_1 == refs_2 && refs_1 == refs_3);
    87.     return 0;
    88. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-18 09:42 | 显示全部楼层
    1. /*
    2. 计算几何
    3. 目录
    4. ㈠ 点的基本运算
    5. 1. 平面上两点之间距离 1
    6. 2. 判断两点是否重合 1
    7. 3. 矢量叉乘 1
    8. 4. 矢量点乘 2
    9. 5. 判断点是否在线段上 2
    10. 6. 求一点饶某点旋转后的坐标 2
    11. 7. 求矢量夹角 2
    12. ㈡ 线段及直线的基本运算
    13. 1. 点与线段的关系 3
    14. 2. 求点到线段所在直线垂线的垂足 4
    15. 3. 点到线段的最近点 4
    16. 4. 点到线段所在直线的距离 4
    17. 5. 点到折线集的最近距离 4
    18. 6. 判断圆是否在多边形内 5
    19. 7. 求矢量夹角余弦 5
    20. 8. 求线段之间的夹角 5
    21. 9. 判断线段是否相交 6
    22. 10.判断线段是否相交但不交在端点处 6
    23. 11.求线段所在直线的方程 6
    24. 12.求直线的斜率 7
    25. 13.求直线的倾斜角 7
    26. 14.求点关于某直线的对称点 7
    27. 15.判断两条直线是否相交及求直线交点 7
    28. 16.判断线段是否相交,如果相交返回交点 7
    29. ㈢ 多边形常用算法模块
    30. 1. 判断多边形是否简单多边形 8
    31. 2. 检查多边形顶点的凸凹性 9
    32. 3. 判断多边形是否凸多边形 9
    33. 4. 求多边形面积 9
    34. 5. 判断多边形顶点的排列方向,方法一 10
    35. 6. 判断多边形顶点的排列方向,方法二 10
    36. 7. 射线法判断点是否在多边形内 10
    37. 8. 判断点是否在凸多边形内 11
    38. 9. 寻找点集的graham算法 12
    39. 10.寻找点集凸包的卷包裹法 13
    40. 11.判断线段是否在多边形内 14
    41. 12.求简单多边形的重心 15
    42. 13.求凸多边形的重心 17
    43. 14.求肯定在给定多边形内的一个点 17
    44. 15.求从多边形外一点出发到该多边形的切线 18
    45. 16.判断多边形的核是否存在 19
    46. ㈣ 圆的基本运算
    47. 1 .点是否在圆内 20
    48. 2 .求不共线的三点所确定的圆 21
    49. ㈤ 矩形的基本运算
    50. 1.已知矩形三点坐标,求第4点坐标 22
    51. ㈥ 常用算法的描述 22
    52. ㈦ 补充
    53. 1.两圆关系: 24
    54. 2.判断圆是否在矩形内: 24
    55. 3.点到平面的距离: 25
    56. 4.点是否在直线同侧: 25
    57. 5.镜面反射线: 25
    58. 6.矩形包含: 26
    59. 7.两圆交点: 27
    60. 8.两圆公共面积: 28
    61. 9. 圆和直线关系: 29
    62. 10. 内切圆: 30
    63. 11. 求切点: 31
    64. 12. 线段的左右旋: 31
    65. 13.公式: 32
    66. */
    67. /* 需要包含的头文件 */
    68. #include <cmath >

    69. /* 常用的常量定义 */
    70. const double        INF                = 1E200   
    71. const double        EP                = 1E-10
    72. const int                MAXV        = 300
    73. const double        PI                = 3.14159265

    74. /* 基本几何结构 */
    75. struct POINT
    76. {
    77.         double x;
    78.         double y;
    79.         POINT(double a=0, double b=0) { x=a; y=b;} //constructor
    80. };
    81. struct LINESEG
    82. {
    83.         POINT s;
    84.         POINT e;
    85.         LINESEG(POINT a, POINT b) { s=a; e=b;}
    86.         LINESEG() { }
    87. };
    88. struct LINE           // 直线的解析方程 a*x+b*y+c=0  为统一表示,约定 a >= 0
    89. {
    90.    double a;
    91.    double b;
    92.    double c;
    93.    LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;}
    94. };

    95. /**********************
    96. *                    *
    97. *   点的基本运算     *
    98. *                    *
    99. **********************/

    100. double dist(POINT p1,POINT p2)                // 返回两点之间欧氏距离
    101. {
    102.         return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );
    103. }
    104. bool equal_point(POINT p1,POINT p2)           // 判断两个点是否重合  
    105. {
    106.         return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) );
    107. }
    108. /******************************************************************************
    109. r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉积
    110. r>0:ep在矢量opsp的逆时针方向;
    111. r=0:opspep三点共线;
    112. r<0:ep在矢量opsp的顺时针方向
    113. *******************************************************************************/
    114. double multiply(POINT sp,POINT ep,POINT op)
    115. {
    116.         return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));
    117. }
    118. /*
    119. r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积,如果两个矢量都非零矢量
    120. r<0:两矢量夹角为锐角;
    121. r=0:两矢量夹角为直角;
    122. r>0:两矢量夹角为钝角
    123. *******************************************************************************/
    124. double dotmultiply(POINT p1,POINT p2,POINT p0)
    125. {
    126.         return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y));
    127. }
    128. /******************************************************************************
    129. 判断点p是否在线段l上
    130. 条件:(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)
    131. *******************************************************************************/
    132. bool online(LINESEG l,POINT p)
    133. {
    134.         return( (multiply(l.e,p,l.s)==0) &&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) );
    135. }
    136. // 返回点p以点o为圆心逆时针旋转alpha(单位:弧度)后所在的位置
    137. POINT rotate(POINT o,double alpha,POINT p)
    138. {
    139.         POINT tp;
    140.         p.x-=o.x;
    141.         p.y-=o.y;
    142.         tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;
    143.         tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;
    144.         return tp;
    145. }
    146. /* 返回顶角在o点,起始边为os,终止边为oe的夹角(单位:弧度)
    147.         角度小于pi,返回正值
    148.         角度大于pi,返回负值
    149.         可以用于求线段之间的夹角
    150. 原理:
    151.         r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))
    152.         r'= multiply(s,e,o)
    153.         r >= 1        angle = 0;
    154.         r <= -1        angle = -PI
    155.         -1<r<1 && r'>0        angle = arccos(r)
    156.         -1<r<1 && r'<=0        angle = -arccos(r)
    157. */
    158. double angle(POINT o,POINT s,POINT e)
    159. {
    160.         double cosfi,fi,norm;
    161.         double dsx = s.x - o.x;
    162.         double dsy = s.y - o.y;
    163.         double dex = e.x - o.x;
    164.         double dey = e.y - o.y;

    165.         cosfi=dsx*dex+dsy*dey;
    166.         norm=(dsx*dsx+dsy*dsy)*(dex*dex+dey*dey);
    167.         cosfi /= sqrt( norm );

    168.         if (cosfi >=  1.0 ) return 0;
    169.         if (cosfi <= -1.0 ) return -3.1415926;

    170.         fi=acos(cosfi);
    171.         if (dsx*dey-dsy*dex>0) return fi;      // 说明矢量os 在矢量 oe的顺时针方向
    172.         return -fi;
    173. }
    174.   /*****************************\
    175.   *                             *
    176.   *      线段及直线的基本运算   *
    177.   *                             *
    178.   \*****************************/

    179. /* 判断点与线段的关系,用途很广泛
    180. 本函数是根据下面的公式写的,P是点C到线段AB所在直线的垂足
    181.                 AC dot AB
    182.         r =     ---------
    183.                  ||AB||^2
    184.              (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    185.           = -------------------------------
    186.                           L^2
    187.     r has the following meaning:
    188.         r=0      P = A
    189.         r=1      P = B
    190.         r<0                 P is on the backward extension of AB
    191.                 r>1      P is on the forward extension of AB
    192.         0<r<1         P is interior to AB
    193. */
    194. double relation(POINT p,LINESEG l)
    195. {
    196.         LINESEG tl;
    197.         tl.s=l.s;
    198.         tl.e=p;
    199.         return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e));
    200. }
    201. // 求点C到线段AB所在直线的垂足 P
    202. POINT perpendicular(POINT p,LINESEG l)
    203. {
    204.         double r=relation(p,l);
    205.         POINT tp;
    206.         tp.x=l.s.x+r*(l.e.x-l.s.x);
    207.         tp.y=l.s.y+r*(l.e.y-l.s.y);
    208.         return tp;
    209. }
    210. /* 求点p到线段l的最短距离,并返回线段上距该点最近的点np
    211. 注意:np是线段l上到点p最近的点,不一定是垂足 */
    212. double ptolinesegdist(POINT p,LINESEG l,POINT &np)
    213. {
    214.         double r=relation(p,l);
    215.         if(r<0)
    216.         {
    217.                 np=l.s;
    218.                 return dist(p,l.s);
    219.         }
    220.         if(r>1)
    221.         {
    222.                 np=l.e;
    223.                 return dist(p,l.e);
    224.         }
    225.         np=perpendicular(p,l);
    226.         return dist(p,np);
    227. }
    228. // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别  
    229. double ptoldist(POINT p,LINESEG l)
    230. {
    231.         return abs(multiply(p,l.e,l.s))/dist(l.s,l.e);
    232. }
    233. /* 计算点到折线集的最近距离,并返回最近点.
    234. 注意:调用的是ptolineseg()函数 */
    235. double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q)
    236. {
    237.         int i;
    238.         double cd=double(INF),td;
    239.         LINESEG l;
    240.         POINT tq,cq;

    241.         for(i=0;i<vcount-1;i++)
    242.         {
    243.                 l.s=pointset[i];

    244.                 l.e=pointset[i+1];
    245.                 td=ptolinesegdist(p,l,tq);
    246.                 if(td<cd)
    247.                 {
    248.                         cd=td;
    249.                         cq=tq;
    250.                 }
    251.         }
    252.         q=cq;
    253.         return cd;
    254. }
    255. /* 判断圆是否在多边形内.ptolineseg()函数的应用2 */
    256. bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[])
    257. {
    258.         POINT q;
    259.         double d;
    260.         q.x=0;
    261.         q.y=0;
    262.         d=ptopointset(vcount,polygon,center,q);
    263.         if(d<radius||fabs(d-radius)<EP)
    264.                 return true;
    265.         else
    266.                 return false;
    267. }
    268. /* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意:如果想从余弦求夹角的话,注意反余弦函数的定义域是从 0到pi */
    269. double cosine(LINESEG l1,LINESEG l2)
    270. {
    271.         return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) +
    272.         (l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) );
    273. }
    274. // 返回线段l1与l2之间的夹角 单位:弧度 范围(-pi,pi)
    275. double lsangle(LINESEG l1,LINESEG l2)
    276. {
    277.         POINT o,s,e;
    278.         o.x=o.y=0;
    279.         s.x=l1.e.x-l1.s.x;
    280.         s.y=l1.e.y-l1.s.y;
    281.         e.x=l2.e.x-l2.s.x;
    282.         e.y=l2.e.y-l2.s.y;
    283.         return angle(o,s,e);
    284. }
    285. // 如果线段u和v相交(包括相交在端点处)时,返回true
    286. //
    287. //判断P1P2跨立Q1Q2的依据是:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
    288. //判断Q1Q2跨立P1P2的依据是:( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。
    289. bool intersect(LINESEG u,LINESEG v)
    290. {
    291.         return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&&                     //排斥实验
    292.                         (max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&&
    293.                         (max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&&
    294.                         (max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&&
    295.                         (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&&         //跨立实验
    296.                         (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0));
    297. }
    298. //  (线段u和v相交)&&(交点不是双方的端点) 时返回true   
    299. bool intersect_A(LINESEG u,LINESEG v)
    300. {
    301.         return        ((intersect(u,v))&&
    302.                         (!online(u,v.s))&&
    303.                         (!online(u,v.e))&&
    304.                         (!online(v,u.e))&&
    305.                         (!online(v,u.s)));
    306. }
    307. // 线段v所在直线与线段u相交时返回true;方法:判断线段u是否跨立线段v  
    308. bool intersect_l(LINESEG u,LINESEG v)
    309. {
    310.         return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0;
    311. }
    312. // 根据已知两点坐标,求过这两点的直线解析方程: a*x+b*y+c = 0  (a >= 0)  
    313. LINE makeline(POINT p1,POINT p2)
    314. {
    315.         LINE tl;
    316.         int sign = 1;
    317.         tl.a=p2.y-p1.y;
    318.         if(tl.a<0)
    319.         {
    320.                 sign = -1;
    321.                 tl.a=sign*tl.a;
    322.         }
    323.         tl.b=sign*(p1.x-p2.x);
    324.         tl.c=sign*(p1.y*p2.x-p1.x*p2.y);
    325.         return tl;
    326. }
    327. // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200
    328. double slope(LINE l)
    329. {
    330.         if(abs(l.a) < 1e-20)
    331.                 return 0;
    332.         if(abs(l.b) < 1e-20)
    333.                 return INF;
    334.         return -(l.a/l.b);
    335. }
    336. // 返回直线的倾斜角alpha ( 0 - pi)
    337. double alpha(LINE l)
    338. {
    339.         if(abs(l.a)< EP)
    340.                 return 0;
    341.         if(abs(l.b)< EP)
    342.                 return PI/2;
    343.         double k=slope(l);
    344.         if(k>0)
    345.                 return atan(k);
    346.         else
    347.                 return PI+atan(k);
    348. }
    349. // 求点p关于直线l的对称点  
    350. POINT symmetry(LINE l,POINT p)
    351. {
    352.    POINT tp;
    353.    tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
    354.    tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
    355.    return tp;
    356. }
    357. // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交,返回true,且返回交点p  
    358. bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1,L2
    359. {
    360.         double d=l1.a*l2.b-l2.a*l1.b;
    361.         if(abs(d)<EP) // 不相交
    362.                 return false;
    363.         p.x = (l2.c*l1.b-l1.c*l2.b)/d;
    364.         p.y = (l2.a*l1.c-l1.a*l2.c)/d;
    365.         return true;
    366. }
    367. // 如果线段l1和l2相交,返回true且交点由(inter)返回,否则返回false
    368. bool intersection(LINESEG l1,LINESEG l2,POINT &inter)
    369. {
    370.         LINE ll1,ll2;
    371.         ll1=makeline(l1.s,l1.e);
    372.         ll2=makeline(l2.s,l2.e);
    373.         if(lineintersect(ll1,ll2,inter))
    374.                 return online(l1,inter);
    375.         else
    376.                 return false;
    377. }

    378. /******************************\
    379. *                                                          *
    380. * 多边形常用算法模块                  *
    381. *                                                          *
    382. \******************************/

    383. // 如果无特别说明,输入多边形顶点要求按逆时针排列

    384. /*
    385. 返回值:输入的多边形是简单多边形,返回true
    386. 要 求:输入顶点序列按逆时针排序
    387. 说 明:简单多边形定义:
    388. 1:循环排序中相邻线段对的交是他们之间共有的单个点
    389. 2:不相邻的线段不相交
    390. 本程序默认第一个条件已经满足
    391. */
    392. bool issimple(int vcount,POINT polygon[])
    393. {
    394.         int i,cn;
    395.         LINESEG l1,l2;
    396.         for(i=0;i<vcount;i++)
    397.         {
    398.                 l1.s=polygon[i];
    399.                 l1.e=polygon[(i+1)%vcount];
    400.                 cn=vcount-3;
    401.                 while(cn)
    402.                 {
    403.                         l2.s=polygon[(i+2)%vcount];
    404.                         l2.e=polygon[(i+3)%vcount];
    405.                         if(intersect(l1,l2))
    406.                                 break;
    407.                         cn--;
    408.                 }
    409.                 if(cn)
    410.                         return false;
    411.         }
    412.         return true;
    413. }
    414. // 返回值:按输入顺序返回多边形顶点的凸凹性判断,bc[i]=1,iff:第i个顶点是凸顶点
    415. void checkconvex(int vcount,POINT polygon[],bool bc[])
    416. {
    417.         int i,index=0;
    418.         POINT tp=polygon[0];
    419.         for(i=1;i<vcount;i++) // 寻找第一个凸顶点
    420.         {
    421.                 if(polygon[i].y<tp.y||(polygon[i].y == tp.y&&polygon[i].x<tp.x))
    422.                 {
    423.                         tp=polygon[i];
    424.                         index=i;
    425.                 }
    426.         }
    427.         int count=vcount-1;
    428.         bc[index]=1;
    429.         while(count) // 判断凸凹性
    430.         {
    431.                 if(multiply(polygon[(index+1)%vcount],polygon[(index+2)%vcount],polygon[index])>=0 )
    432.                         bc[(index+1)%vcount]=1;
    433.                 else
    434.                         bc[(index+1)%vcount]=0;
    435.                 index++;
    436.                 count--;
    437.         }
    438. }
    439. // 返回值:多边形polygon是凸多边形时,返回true  
    440. bool isconvex(int vcount,POINT polygon[])
    441. {
    442.         bool bc[MAXV];
    443.         checkconvex(vcount,polygon,bc);
    444.         for(int i=0;i<vcount;i++) // 逐一检查顶点,是否全部是凸顶点
    445.                 if(!bc[i])
    446.                         return false;
    447.         return true;
    448. }
    449. // 返回多边形面积(signed);输入顶点按逆时针排列时,返回正值;否则返回负值
    450. double area_of_polygon(int vcount,POINT polygon[])
    451. {
    452.         int i;
    453.         double s;
    454.         if (vcount<3)
    455.                 return 0;
    456.         s=polygon[0].y*(polygon[vcount-1].x-polygon[1].x);
    457.         for (i=1;i<vcount;i++)
    458.                 s+=polygon[i].y*(polygon[(i-1)].x-polygon[(i+1)%vcount].x);
    459.         return s/2;
    460. }
    461. // 如果输入顶点按逆时针排列,返回true
    462. bool isconterclock(int vcount,POINT polygon[])
    463. {
    464.         return area_of_polygon(vcount,polygon)>0;
    465. }
    466. // 另一种判断多边形顶点排列方向的方法  
    467. bool isccwize(int vcount,POINT polygon[])
    468. {
    469.         int i,index;
    470.         POINT a,b,v;
    471.         v=polygon[0];
    472.         index=0;
    473.         for(i=1;i<vcount;i++) // 找到最低且最左顶点,肯定是凸顶点
    474.         {
    475.                 if(polygon[i].y<v.y||polygon[i].y == v.y && polygon[i].x<v.x)
    476.                 {
    477.                         index=i;
    478.                 }
    479.         }
    480.         a=polygon[(index-1+vcount)%vcount]; // 顶点v的前一顶点
    481.         b=polygon[(index+1)%vcount]; // 顶点v的后一顶点
    482.         return multiply(v,b,a)>0;
    483. }
    484. /********************************************************************************************   
    485. 射线法判断点q与多边形polygon的位置关系,要求polygon为简单多边形,顶点逆时针排列
    486.    如果点在多边形内:   返回0
    487.    如果点在多边形边上: 返回1
    488.    如果点在多边形外:        返回2
    489. *********************************************************************************************/
    490. int insidepolygon(int vcount,POINT Polygon[],POINT q)
    491. {
    492.         int c=0,i,n;
    493.         LINESEG l1,l2;
    494.         bool bintersect_a,bonline1,bonline2,bonline3;
    495.         double r1,r2;

    496.         l1.s=q;
    497.         l1.e=q;
    498.         l1.e.x=double(INF);
    499.         n=vcount;
    500.         for (i=0;i<vcount;i++)
    501.         {
    502.                 l2.s=Polygon[i];
    503.                 l2.e=Polygon[(i+1)%n];
    504.                 if(online(l2,q))
    505.                         return 1; // 如果点在边上,返回1
    506.                 if ( (bintersect_a=intersect_A(l1,l2))|| // 相交且不在端点
    507.                 ( (bonline1=online(l1,Polygon[(i+1)%n]))&& // 第二个端点在射线上
    508.                 ( (!(bonline2=online(l1,Polygon[(i+2)%n])))&& /* 前一个端点和后一个端点在射线两侧 */
    509.                 ((r1=multiply(Polygon[i],Polygon[(i+1)%n],l1.s)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.s))>0) ||   
    510.                 (bonline3=online(l1,Polygon[(i+2)%n]))&&     /* 下一条边是水平线,前一个端点和后一个端点在射线两侧  */
    511.                         ((r2=multiply(Polygon[i],Polygon[(i+2)%n],l1.s)*multiply(Polygon[(i+2)%n],
    512.                 Polygon[(i+3)%n],l1.s))>0)
    513.                                 )
    514.                         )
    515.                 ) c++;
    516.         }
    517.         if(c%2 == 1)
    518.                 return 0;
    519.         else
    520.                 return 2;
    521. }
    522. //点q是凸多边形polygon内时,返回true;注意:多边形polygon一定要是凸多边形  
    523. bool InsideConvexPolygon(int vcount,POINT polygon[],POINT q) // 可用于三角形!
    524. {
    525.         POINT p;
    526.         LINESEG l;
    527.         int i;
    528.         p.x=0;p.y=0;
    529.         for(i=0;i<vcount;i++) // 寻找一个肯定在多边形polygon内的点p:多边形顶点平均值
    530.         {
    531.                 p.x+=polygon[i].x;
    532.                 p.y+=polygon[i].y;
    533.         }
    534.         p.x /= vcount;
    535.         p.y /= vcount;

    536.         for(i=0;i<vcount;i++)
    537.         {
    538.                 l.s=polygon[i];l.e=polygon[(i+1)%vcount];
    539.                 if(multiply(p,l.e,l.s)*multiply(q,l.e,l.s)<0) /* 点p和点q在边l的两侧,说明点q肯定在多边形外 */
    540.                 break;
    541.         }
    542.         return (i==vcount);
    543. }
    544. /**********************************************
    545. 寻找凸包的graham 扫描法
    546. PointSet为输入的点集;
    547. ch为输出的凸包上的点集,按照逆时针方向排列;
    548. n为PointSet中的点的数目
    549. len为输出的凸包上的点的个数
    550. **********************************************/
    551. void Graham_scan(POINT PointSet[],POINT ch[],int n,int &len)
    552. {
    553.         int i,j,k=0,top=2;
    554.         POINT tmp;
    555.         // 选取PointSet中y坐标最小的点PointSet[k],如果这样的点有多个,则取最左边的一个
    556.         for(i=1;i<n;i++)
    557.                 if ( PointSet[i].y<PointSet[k].y || (PointSet[i].y==PointSet[k].y) && (PointSet[i].x<PointSet[k].x) )
    558.                         k=i;
    559.         tmp=PointSet[0];
    560.         PointSet[0]=PointSet[k];
    561.         PointSet[k]=tmp; // 现在PointSet中y坐标最小的点在PointSet[0]
    562.         for (i=1;i<n-1;i++) /* 对顶点按照相对PointSet[0]的极角从小到大进行排序,极角相同的按照距离PointSet[0]从近到远进行排序 */
    563.         {
    564.                 k=i;
    565.                 for (j=i+1;j<n;j++)
    566.                         if ( multiply(PointSet[j],PointSet[k],PointSet[0])>0 ||  // 极角更小   
    567.                                 (multiply(PointSet[j],PointSet[k],PointSet[0])==0) && /* 极角相等,距离更短 */        
    568.                                 dist(PointSet[0],PointSet[j])<dist(PointSet[0],PointSet[k])
    569.                            )
    570.                                 k=j;
    571.                 tmp=PointSet[i];
    572.                 PointSet[i]=PointSet[k];
    573.                 PointSet[k]=tmp;
    574.         }
    575.         ch[0]=PointSet[0];
    576.         ch[1]=PointSet[1];
    577.         ch[2]=PointSet[2];
    578.         for (i=3;i<n;i++)
    579.         {
    580.                 while (multiply(PointSet[i],ch[top],ch[top-1])>=0)
    581.                         top--;
    582.                 ch[++top]=PointSet[i];
    583.         }
    584.         len=top+1;
    585. }
    586. // 卷包裹法求点集凸壳,参数说明同graham算法   
    587. void ConvexClosure(POINT PointSet[],POINT ch[],int n,int &len)
    588. {
    589.         int top=0,i,index,first;
    590.         double curmax,curcos,curdis;
    591.         POINT tmp;
    592.         LINESEG l1,l2;
    593.         bool use[MAXV];
    594.         tmp=PointSet[0];
    595.         index=0;
    596.         // 选取y最小点,如果多于一个,则选取最左点
    597.         for(i=1;i<n;i++)
    598.         {
    599.                 if(PointSet[i].y<tmp.y||PointSet[i].y == tmp.y&&PointSet[i].x<tmp.x)
    600.                 {
    601.                         index=i;
    602.                 }
    603.                 use[i]=false;
    604.         }
    605.         tmp=PointSet[index];
    606.         first=index;
    607.         use[index]=true;

    608.         index=-1;
    609.         ch[top++]=tmp;
    610.         tmp.x-=100;
    611.         l1.s=tmp;
    612.         l1.e=ch[0];
    613.         l2.s=ch[0];

    614.         while(index!=first)
    615.         {
    616.                 curmax=-100;
    617.                 curdis=0;
    618.                 // 选取与最后一条确定边夹角最小的点,即余弦值最大者
    619.                 for(i=0;i<n;i++)
    620.                 {
    621.                         if(use[i])continue;
    622.                         l2.e=PointSet[i];
    623.                         curcos=cosine(l1,l2); // 根据cos值求夹角余弦,范围在 (-1 -- 1 )
    624.                         if(curcos>curmax || fabs(curcos-curmax)<1e-6 && dist(l2.s,l2.e)>curdis)
    625.                         {
    626.                                 curmax=curcos;
    627.                                 index=i;
    628.                                 curdis=dist(l2.s,l2.e);
    629.                         }
    630.                 }
    631.                 use[first]=false;            //清空第first个顶点标志,使最后能形成封闭的hull
    632.                 use[index]=true;
    633.                 ch[top++]=PointSet[index];
    634.                 l1.s=ch[top-2];
    635.                 l1.e=ch[top-1];
    636.                 l2.s=ch[top-1];
    637.         }
    638.         len=top-1;
    639. }
    640. /*********************************************************************************************  
    641.         判断线段是否在简单多边形内(注意:如果多边形是凸多边形,下面的算法可以化简)
    642.            必要条件一:线段的两个端点都在多边形内;
    643.         必要条件二:线段和多边形的所有边都不内交;
    644.         用途:        1. 判断折线是否在简单多边形内
    645.                         2. 判断简单多边形是否在另一个简单多边形内
    646. **********************************************************************************************/
    647. bool LinesegInsidePolygon(int vcount,POINT polygon[],LINESEG l)
    648. {
    649.         // 判断线端l的端点是否不都在多边形内
    650.         if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e))
    651.                 return false;
    652.         int top=0,i,j;
    653.         POINT PointSet[MAXV],tmp;
    654.         LINESEG s;

    655.         for(i=0;i<vcount;i++)
    656.         {
    657.                 s.s=polygon[i];
    658.                 s.e=polygon[(i+1)%vcount];
    659.                 if(online(s,l.s)) //线段l的起始端点在线段s上
    660.                         PointSet[top++]=l.s;
    661.                 else if(online(s,l.e)) //线段l的终止端点在线段s上
    662.                         PointSet[top++]=l.e;
    663.                 else
    664.                 {
    665.                         if(online(l,s.s)) //线段s的起始端点在线段l上
    666.                                 PointSet[top++]=s.s;
    667.                         else if(online(l,s.e)) // 线段s的终止端点在线段l上
    668.                                 PointSet[top++]=s.e;
    669.                         else
    670.                         {
    671.                                 if(intersect(l,s)) // 这个时候如果相交,肯定是内交,返回false
    672.                                 return false;
    673.                         }
    674.                 }
    675.         }

    676.         for(i=0;i<top-1;i++) /* 冒泡排序,x坐标小的排在前面;x坐标相同者,y坐标小的排在前面 */
    677.         {
    678.                 for(j=i+1;j<top;j++)
    679.                 {
    680.                         if( PointSet[i].x>PointSet[j].x || fabs(PointSet[i].x-PointSet[j].x)<EP && PointSet[i].y>PointSet[j].y )
    681.                         {
    682.                                 tmp=PointSet[i];
    683.                                 PointSet[i]=PointSet[j];
    684.                                 PointSet[j]=tmp;
    685.                         }
    686.                 }
    687.         }

    688.         for(i=0;i<top-1;i++)
    689.         {
    690.                 tmp.x=(PointSet[i].x+PointSet[i+1].x)/2; //得到两个相邻交点的中点
    691.                 tmp.y=(PointSet[i].y+PointSet[i+1].y)/2;
    692.                 if(!insidepolygon(vcount,polygon,tmp))
    693.                         return false;
    694.         }
    695.         return true;
    696. }
    697. /*********************************************************************************************  
    698. 求任意简单多边形polygon的重心
    699. 需要调用下面几个函数:
    700.         void AddPosPart(); 增加右边区域的面积
    701.         void AddNegPart(); 增加左边区域的面积
    702.         void AddRegion(); 增加区域面积
    703. 在使用该程序时,如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简,
    704. 但要注意函数的声明和调用要做相应变化
    705. **********************************************************************************************/
    706. void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr)
    707. {
    708.         if (abs(wtr + w)<1e-10 ) return; // detect zero regions
    709.         xtr = ( wtr*xtr + w*x ) / ( wtr + w );
    710.         ytr = ( wtr*ytr + w*y ) / ( wtr + w );
    711.         wtr = w + wtr;
    712.         return;
    713. }
    714. void AddNegPart(double x, ouble y, double w, double &xtl, double &ytl, double &wtl)
    715. {
    716.         if ( abs(wtl + w)<1e-10 )
    717.                 return; // detect zero regions

    718.         xtl = ( wtl*xtl + w*x ) / ( wtl + w );
    719.         ytl = ( wtl*ytl + w*y ) / ( wtl + w );
    720.         wtl = w + wtl;
    721.         return;
    722. }
    723. void AddRegion ( double x1, double y1, double x2, double y2, double &xtr, double &ytr,
    724.                 double &wtr, double &xtl, double &ytl, double &wtl )
    725. {
    726.         if ( abs (x1 - x2)< 1e-10 )
    727.                 return;

    728.         if ( x2 > x1 )
    729.         {
    730.                 AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全局变量变化处 */
    731.                 AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr);   
    732.                 // triangle 全局变量变化处
    733.         }
    734.         else
    735.         {
    736.                 AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl);  
    737.                 // rectangle 全局变量变化处
    738.                 AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl);  
    739.                 // triangle  全局变量变化处
    740.         }
    741. }
    742. POINT cg_simple(int vcount,POINT polygon[])
    743. {
    744.         double xtr,ytr,wtr,xtl,ytl,wtl;        
    745.         //注意: 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去
    746.         POINT p1,p2,tp;
    747.         xtr = ytr = wtr = 0.0;
    748.         xtl = ytl = wtl = 0.0;
    749.         for(int i=0;i<vcount;i++)
    750.         {
    751.                 p1=polygon[i];
    752.                 p2=polygon[(i+1)%vcount];
    753.                 AddRegion(p1.x,p1.y,p2.x,p2.y,xtr,ytr,wtr,xtl,ytl,wtl); //全局变量变化处
    754.         }
    755.         tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl);
    756.         tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl);
    757.         return tp;
    758. }
    759. // 求凸多边形的重心,要求输入多边形按逆时针排序
    760. POINT gravitycenter(int vcount,POINT polygon[])
    761. {
    762.         POINT tp;
    763.         double x,y,s,x0,y0,cs,k;
    764.         x=0;y=0;s=0;
    765.         for(int i=1;i<vcount-1;i++)
    766.         {
    767.                 x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3;
    768.                 y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求当前三角形的重心
    769.                 cs=multiply(polygon[i],polygon[i+1],polygon[0])/2;
    770.                 //三角形面积可以直接利用该公式求解
    771.                 if(abs(s)<1e-20)
    772.                 {
    773.                         x=x0;y=y0;s+=cs;continue;
    774.                 }
    775.                 k=cs/s; //求面积比例
    776.                 x=(x+k*x0)/(1+k);
    777.                 y=(y+k*y0)/(1+k);
    778.                 s += cs;
    779.         }
    780.         tp.x=x;
    781.         tp.y=y;
    782.         return tp;
    783. }

    784. /************************************************
    785. 给定一简单多边形,找出一个肯定在该多边形内的点
    786. 定理1        :每个多边形至少有一个凸顶点
    787. 定理2        :顶点数>=4的简单多边形至少有一条对角线
    788. 结论        : x坐标最大,最小的点肯定是凸顶点
    789.         y坐标最大,最小的点肯定是凸顶点            
    790. ************************************************/
    791. POINT a_point_insidepoly(int vcount,POINT polygon[])
    792. {
    793.         POINT v,a,b,r;
    794.         int i,index;
    795.         v=polygon[0];
    796.         index=0;
    797.         for(i=1;i<vcount;i++) //寻找一个凸顶点
    798.         {
    799.                 if(polygon[i].y<v.y)
    800.                 {
    801.                         v=polygon[i];
    802.                         index=i;
    803.                 }
    804.         }
    805.         a=polygon[(index-1+vcount)%vcount]; //得到v的前一个顶点
    806.         b=polygon[(index+1)%vcount]; //得到v的后一个顶点
    807.         POINT tri[3],q;
    808.         tri[0]=a;tri[1]=v;tri[2]=b;
    809.         double md=INF;
    810.         int in1=index;
    811.         bool bin=false;
    812.         for(i=0;i<vcount;i++) //寻找在三角形avb内且离顶点v最近的顶点q
    813.         {
    814.                 if(i == index)continue;
    815.                 if(i == (index-1+vcount)%vcount)continue;
    816.                 if(i == (index+1)%vcount)continue;
    817.                 if(!InsideConvexPolygon(3,tri,polygon[i]))continue;
    818.                 bin=true;
    819.                 if(dist(v,polygon[i])<md)
    820.                 {
    821.                         q=polygon[i];
    822.                         md=dist(v,q);
    823.                 }
    824.         }
    825.         if(!bin) //没有顶点在三角形avb内,返回线段ab中点
    826.         {
    827.                 r.x=(a.x+b.x)/2;
    828.                 r.y=(a.y+b.y)/2;
    829.                 return r;
    830.         }
    831.         r.x=(v.x+q.x)/2; //返回线段vq的中点
    832.         r.y=(v.y+q.y)/2;
    833.         return r;
    834. }
    835. /***********************************************************************************************
    836. 求从多边形外一点p出发到一个简单多边形的切线,如果存在返回切点,其中rp点是右切点,lp是左切点
    837. 注意:p点一定要在多边形外 ,输入顶点序列是逆时针排列
    838. 原 理:        如果点在多边形内肯定无切线;凸多边形有唯一的两个切点,凹多边形就可能有多于两个的切点;
    839.                 如果polygon是凸多边形,切点只有两个只要找到就可以,可以化简此算法
    840.                 如果是凹多边形还有一种算法可以求解:先求凹多边形的凸包,然后求凸包的切线
    841. /***********************************************************************************************/
    842. void pointtangentpoly(int vcount,POINT polygon[],POINT p,POINT &rp,POINT &lp)
    843. {
    844.         LINESEG ep,en;
    845.         bool blp,bln;
    846.         rp=polygon[0];
    847.         lp=polygon[0];
    848.         for(int i=1;i<vcount;i++)
    849.         {
    850.                 ep.s=polygon[(i+vcount-1)%vcount];
    851.                 ep.e=polygon[i];
    852.                 en.s=polygon[i];
    853.                 en.e=polygon[(i+1)%vcount];
    854.                 blp=multiply(ep.e,p,ep.s)>=0;                // p is to the left of pre edge
    855.                 bln=multiply(en.e,p,en.s)>=0;                // p is to the left of next edge
    856.                 if(!blp&&bln)
    857.                 {
    858.                         if(multiply(polygon[i],rp,p)>0)           // polygon[i] is above rp
    859.                         rp=polygon[i];
    860.                 }
    861.                 if(blp&&!bln)
    862.                 {
    863.                         if(multiply(lp,polygon[i],p)>0)           // polygon[i] is below lp
    864.                         lp=polygon[i];
    865.                 }
    866.         }
    867.         return ;
    868. }
    869. // 如果多边形polygon的核存在,返回true,返回核上的一点p.顶点按逆时针方向输入  
    870. bool core_exist(int vcount,POINT polygon[],POINT &p)
    871. {
    872.         int i,j,k;
    873.         LINESEG l;
    874.         LINE lineset[MAXV];
    875.         for(i=0;i<vcount;i++)
    876.         {
    877.                 lineset[i]=makeline(polygon[i],polygon[(i+1)%vcount]);
    878.         }
    879.         for(i=0;i<vcount;i++)
    880.         {
    881.                 for(j=0;j<vcount;j++)
    882.                 {
    883.                         if(i == j)continue;
    884.                         if(lineintersect(lineset[i],lineset[j],p))
    885.                         {
    886.                                 for(k=0;k<vcount;k++)
    887.                                 {
    888.                                         l.s=polygon[k];
    889.                                         l.e=polygon[(k+1)%vcount];
    890.                                         if(multiply(p,l.e,l.s)>0)      
    891.                                         //多边形顶点按逆时针方向排列,核肯定在每条边的左侧或边上
    892.                                         break;
    893.                                 }
    894.                                 if(k == vcount)             //找到了一个核上的点
    895.                                 break;
    896.                         }
    897.                 }
    898.                 if(j<vcount) break;
    899.         }
    900.         if(i<vcount)
    901.                 return true;
    902.         else
    903.                 return false;
    904. }
    905. /*************************\
    906. *                                                 *
    907. * 圆的基本运算           *
    908. *                                             *
    909. \*************************/
    910. /******************************************************************************
    911. 返回值        : 点p在圆内(包括边界)时,返回true
    912. 用途        : 因为圆为凸集,所以判断点集,折线,多边形是否在圆内时,
    913.         只需要逐一判断点是否在圆内即可。
    914. *******************************************************************************/
    915. bool point_in_circle(POINT o,double r,POINT p)
    916. {
    917.         double d2=(p.x-o.x)*(p.x-o.x)+(p.y-o.y)*(p.y-o.y);
    918.         double r2=r*r;
    919.         return d2<r2||abs(d2-r2)<EP;
    920. }
    921. /******************************************************************************
    922. 用 途        :求不共线的三点确定一个圆
    923. 输 入        :三个点p1,p2,p3
    924. 返回值        :如果三点共线,返回false;反之,返回true。圆心由q返回,半径由r返回
    925. *******************************************************************************/
    926. bool cocircle(POINT p1,POINT p2,POINT p3,POINT &q,double &r)
    927. {
    928.         double x12=p2.x-p1.x;
    929.         double y12=p2.y-p1.y;
    930.         double x13=p3.x-p1.x;
    931.         double y13=p3.y-p1.y;
    932.         double z2=x12*(p1.x+p2.x)+y12*(p1.y+p2.y);
    933.         double z3=x13*(p1.x+p3.x)+y13*(p1.y+p3.y);
    934.         double d=2.0*(x12*(p3.y-p2.y)-y12*(p3.x-p2.x));
    935.         if(abs(d)<EP) //共线,圆不存在
    936.                 return false;
    937.         q.x=(y13*z2-y12*z3)/d;
    938.         q.y=(x12*z3-x13*z2)/d;
    939.         r=dist(p1,q);
    940.         return true;
    941. }
    942. int line_circle(LINE l,POINT o,double r,POINT &p1,POINT &p2)
    943. {
    944.         return true;
    945. }

    946. /**************************\
    947. *                                                  *
    948. * 矩形的基本运算          *
    949. *                         *
    950. \**************************/
    951. /*
    952. 说明:因为矩形的特殊性,常用算法可以化简:
    953. 1.判断矩形是否包含点
    954. 只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。
    955. 2.判断线段、折线、多边形是否在矩形中
    956. 因为矩形是个凸集,所以只要判断所有端点是否都在矩形中就可以了。
    957. 3.判断圆是否在矩形中
    958. 圆在矩形中的充要条件是:圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。
    959. */
    960. // 已知矩形的三个顶点(a,b,c),计算第四个顶点d的坐标. 注意:已知的三个顶点可以是无序的
    961. POINT rect4th(POINT a,POINT b,POINT c)
    962. {
    963.         POINT d;
    964.         if(abs(dotmultiply(a,b,c))<EP) // 说明c点是直角拐角处
    965.         {
    966.                 d.x=a.x+b.x-c.x;
    967.                 d.y=a.y+b.y-c.y;
    968.         }
    969.         if(abs(dotmultiply(a,c,b))<EP) // 说明b点是直角拐角处
    970.         {
    971.                 d.x=a.x+c.x-b.x;
    972.                 d.y=a.y+c.y-b.x;
    973.         }
    974.         if(abs(dotmultiply(c,b,a))<EP) // 说明a点是直角拐角处
    975.         {
    976.                 d.x=c.x+b.x-a.x;
    977.                 d.y=c.y+b.y-a.y;
    978.         }
    979.         return d;
    980. }

    981. /*************************\
    982. *                                                *
    983. * 常用算法的描述                *
    984. *                                                *
    985. \*************************/
    986. /*
    987. 尚未实现的算法:
    988. 1. 求包含点集的最小圆
    989. 2. 求多边形的交
    990. 3. 简单多边形的三角剖分
    991. 4. 寻找包含点集的最小矩形
    992. 5. 折线的化简
    993. 6. 判断矩形是否在矩形中
    994. 7. 判断矩形能否放在矩形中
    995. 8. 矩形并的面积与周长
    996. 9. 矩形并的轮廓
    997. 10.矩形并的闭包
    998. 11.矩形的交
    999. 12.点集中的最近点对
    1000. 13.多边形的并
    1001. 14.圆的交与并
    1002. 15.直线与圆的关系
    1003. 16.线段与圆的关系
    1004. 17.求多边形的核监视摄象机
    1005. 18.求点集中不相交点对 railwai
    1006. *//*
    1007. 寻找包含点集的最小矩形
    1008. 原理:该矩形至少一条边与点集的凸壳的某条边共线
    1009. First take the convex hull of the points. Let the resulting convex
    1010. polygon be P. It has been known for some time that the minimum
    1011. area rectangle enclosing P must have one rectangle side flush with
    1012. (i.e., collinear with and overlapping) one edge of P. This geometric
    1013. fact was used by Godfried Toussaint to develop the "rotating calipers"
    1014. algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems
    1015. with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm
    1016. rotates a surrounding rectangle from one flush edge to the next,
    1017. keeping track of the minimum area for each edge. It achieves O(n)
    1018. time (after hull computation). See the "Rotating Calipers Homepage"
    1019. http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description
    1020. and applet.
    1021. *//*
    1022. 折线的化简 伪码如下:
    1023. Input: tol = the approximation tolerance
    1024. L = {V0,V1,...,Vn-1} is any n-vertex polyline
    1025. Set start = 0;
    1026. Set k = 0;
    1027. Set W0 = V0;
    1028. for each vertex Vi (i=1,n-1)
    1029. {
    1030. if Vi is within tol from Vstart
    1031. then ignore it, and continue with the next vertex
    1032. Vi is further than tol away from Vstart
    1033. so add it as a new vertex of the reduced polyline
    1034. Increment k++;
    1035. Set Wk = Vi;
    1036. Set start = i; as the new initial vertex
    1037. }
    1038. Output: W = {W0,W1,...,Wk-1} = the k-vertex simplified polyline
    1039. */
    1040. /********************\
    1041. *                                    *
    1042. * 补充                                *
    1043. *                                        *
    1044. \********************/

    1045. //两圆关系:
    1046. /* 两圆:
    1047. 相离: return 1;
    1048. 外切: return 2;
    1049. 相交: return 3;
    1050. 内切: return 4;
    1051. 内含: return 5;
    1052. */
    1053. int CircleRelation(POINT p1, double r1, POINT p2, double r2)
    1054. {
    1055.         double d = sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) );

    1056.         if( fabs(d-r1-r2) < EP ) // 必须保证前两个if先被判定!
    1057.                 return 2;
    1058.         if( fabs(d-fabs(r1-r2)) < EP )
    1059.                 return 4;
    1060.         if( d > r1+r2 )
    1061.                 return 1;
    1062.         if( d < fabs(r1-r2) )
    1063.                 return 5;
    1064.         if( fabs(r1-r2) < d && d < r1+r2 )
    1065.                 return 3;
    1066.         return 0; // indicate an error!
    1067. }
    1068. //判断圆是否在矩形内:
    1069. // 判定圆是否在矩形内,是就返回true(设矩形水平,且其四个顶点由左上开始按顺时针排列)
    1070. // 调用ptoldist函数,在第4页
    1071. bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4)
    1072. {
    1073.         if( pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y )
    1074.         {
    1075.                 LINESEG line1(pr1, pr2);
    1076.                 LINESEG line2(pr2, pr3);
    1077.                 LINESEG line3(pr3, pr4);
    1078.                 LINESEG line4(pr4, pr1);
    1079.                 if( r<ptoldist(pc,line1) && r<ptoldist(pc,line2) &&        r<ptoldist(pc,line3) && r<ptoldist(pc,line4) )
    1080.                         return true;
    1081.         }
    1082.         return false;
    1083. }
    1084. //点到平面的距离:
    1085. //点到平面的距离,平面用一般式表示ax+by+cz+d=0
    1086. double P2planeDist(double x, double y, double z, double a, double b, double c, double d)
    1087. {
    1088.         return fabs(a*x+b*y+c*z+d) / sqrt(a*a+b*b+c*c);
    1089. }
    1090. //点是否在直线同侧:
    1091. //两个点是否在直线同侧,是则返回true
    1092. bool SameSide(POINT p1, POINT p2, LINE line)
    1093. {
    1094.         return (line.a * p1.x + line.b * p1.y + line.c) *
    1095.         (line.a * p2.x + line.b * p2.y + line.c) > 0;
    1096. }
    1097. //镜面反射线:
    1098. // 已知入射线、镜面,求反射线。
    1099. // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数;  
    1100. //a2,b2,c2为入射光直线方程系数;  
    1101. //a,b,c为反射光直线方程系数.
    1102. // 光是有方向的,使用时注意:入射光向量:<-b2,a2>;反射光向量:<b,-a>.
    1103. // 不要忘记结果中可能会有"negative zeros"
    1104. void reflect(double a1,double b1,double c1,double a2,double b2,double c2,double &a,double &b,double &c)
    1105. {
    1106.         double n,m;
    1107.         double tpb,tpa;
    1108.         tpb=b1*b2+a1*a2;
    1109.         tpa=a2*b1-a1*b2;
    1110.         m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1);
    1111.         n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1);
    1112.         if(fabs(a1*b2-a2*b1)<1e-20)
    1113.         {
    1114.                 a=a2;b=b2;c=c2;
    1115.                 return;
    1116.         }
    1117.         double xx,yy; //(xx,yy)是入射线与镜面的交点。
    1118.         xx=(b1*c2-b2*c1)/(a1*b2-a2*b1);
    1119.         yy=(a2*c1-a1*c2)/(a1*b2-a2*b1);
    1120.         a=n;
    1121.         b=-m;
    1122.         c=m*yy-xx*n;
    1123. }
    1124. //矩形包含:
    1125. // 矩形2(C,D)是否在1(A,B)内
    1126. bool r2inr1(double A,double B,double C,double D)  
    1127. {
    1128.         double X,Y,L,K,DMax;
    1129.         if (A < B)
    1130.         {
    1131.                 double tmp = A;
    1132.                 A = B;
    1133.                 B = tmp;
    1134.         }
    1135.         if (C < D)
    1136.         {
    1137.                 double tmp = C;
    1138.                 C = D;
    1139.                 D = tmp;
    1140.         }
    1141.         if (A > C && B > D)                 // trivial case  
    1142.                 return true;
    1143.         else
    1144.                 if (D >= B)
    1145.                         return false;
    1146.                 else
    1147.                 {
    1148.                         X = sqrt(A * A + B * B);         // outer rectangle's diagonal  
    1149.                         Y = sqrt(C * C + D * D);         // inner rectangle's diagonal  
    1150.                         if (Y < B) // check for marginal conditions
    1151.                                 return true; // the inner rectangle can freely rotate inside
    1152.                         else
    1153.                                 if (Y > X)
    1154.                                         return false;
    1155.                                 else
    1156.                                 {
    1157.                                         L = (B - sqrt(Y * Y - A * A)) / 2;
    1158.                                         K = (A - sqrt(Y * Y - B * B)) / 2;
    1159.                                         DMax = sqrt(L * L + K * K);
    1160.                                         if (D >= DMax)
    1161.                                         return false;
    1162.                                         else
    1163.                                         return true;
    1164.                                 }
    1165.                 }
    1166. }
    1167. //两圆交点:
    1168. // 两圆已经相交(相切)
    1169. void  c2point(POINT p1,double r1,POINT p2,double r2,POINT &rp1,POINT &rp2)
    1170. {
    1171.         double a,b,r;
    1172.         a=p2.x-p1.x;
    1173.         b=p2.y-p1.y;
    1174.         r=(a*a+b*b+r1*r1-r2*r2)/2;
    1175.         if(a==0&&b!=0)
    1176.         {
    1177.                 rp1.y=rp2.y=r/b;
    1178.                 rp1.x=sqrt(r1*r1-rp1.y*rp1.y);
    1179.                 rp2.x=-rp1.x;
    1180.         }
    1181.         else if(a!=0&&b==0)
    1182.         {
    1183.                 rp1.x=rp2.x=r/a;
    1184.                 rp1.y=sqrt(r1*r1-rp1.x*rp2.x);
    1185.                 rp2.y=-rp1.y;
    1186.         }
    1187.         else if(a!=0&&b!=0)
    1188.         {
    1189.                 double delta;
    1190.                 delta=b*b*r*r-(a*a+b*b)*(r*r-r1*r1*a*a);
    1191.                 rp1.y=(b*r+sqrt(delta))/(a*a+b*b);
    1192.                 rp2.y=(b*r-sqrt(delta))/(a*a+b*b);
    1193.                 rp1.x=(r-b*rp1.y)/a;
    1194.                 rp2.x=(r-b*rp2.y)/a;
    1195.         }

    1196.         rp1.x+=p1.x;
    1197.         rp1.y+=p1.y;
    1198.         rp2.x+=p1.x;
    1199.         rp2.y+=p1.y;
    1200. }
    1201. //两圆公共面积:
    1202. // 必须保证相交
    1203. double c2area(POINT p1,double r1,POINT p2,double r2)
    1204. {
    1205.         POINT rp1,rp2;
    1206.         c2point(p1,r1,p2,r2,rp1,rp2);

    1207.         if(r1>r2) //保证r2>r1
    1208.         {
    1209.                 swap(p1,p2);
    1210.                 swap(r1,r2);
    1211.         }
    1212.         double a,b,rr;
    1213.         a=p1.x-p2.x;
    1214.         b=p1.y-p2.y;
    1215.         rr=sqrt(a*a+b*b);

    1216.         double dx1,dy1,dx2,dy2;
    1217.         double sita1,sita2;
    1218.         dx1=rp1.x-p1.x;
    1219.         dy1=rp1.y-p1.y;
    1220.         dx2=rp2.x-p1.x;
    1221.         dy2=rp2.y-p1.y;
    1222.         sita1=acos((dx1*dx2+dy1*dy2)/r1/r1);

    1223.         dx1=rp1.x-p2.x;
    1224.         dy1=rp1.y-p2.y;
    1225.         dx2=rp2.x-p2.x;
    1226.         dy2=rp2.y-p2.y;
    1227.         sita2=acos((dx1*dx2+dy1*dy2)/r2/r2);
    1228.         double s=0;
    1229.         if(rr<r2)//相交弧为优弧
    1230.                 s=r1*r1*(PI-sita1/2+sin(sita1)/2)+r2*r2*(sita2-sin(sita2))/2;
    1231.         else//相交弧为劣弧
    1232.                 s=(r1*r1*(sita1-sin(sita1))+r2*r2*(sita2-sin(sita2)))/2;

    1233.         return s;
    1234. }
    1235. //圆和直线关系:
    1236. //0----相离 1----相切 2----相交
    1237. int clpoint(POINT p,double r,double a,double b,double c,POINT &rp1,POINT &rp2)
    1238. {
    1239.         int res=0;

    1240.         c=c+a*p.x+b*p.y;
    1241.         double tmp;
    1242.         if(a==0&&b!=0)
    1243.         {
    1244.                 tmp=-c/b;
    1245.                 if(r*r<tmp*tmp)
    1246.                         res=0;
    1247.                 else if(r*r==tmp*tmp)
    1248.                 {
    1249.                         res=1;
    1250.                         rp1.y=tmp;
    1251.                         rp1.x=0;
    1252.                 }
    1253.                 else
    1254.                 {
    1255.                         res=2;
    1256.                         rp1.y=rp2.y=tmp;
    1257.                         rp1.x=sqrt(r*r-tmp*tmp);
    1258.                         rp2.x=-rp1.x;
    1259.                 }
    1260.         }
    1261.         else if(a!=0&&b==0)
    1262.         {
    1263.                 tmp=-c/a;
    1264.                 if(r*r<tmp*tmp)
    1265.                         res=0;
    1266.                 else if(r*r==tmp*tmp)
    1267.                 {
    1268.                         res=1;
    1269.                         rp1.x=tmp;
    1270.                         rp1.y=0;
    1271.                 }
    1272.                 else
    1273.                 {
    1274.                         res=2;
    1275.                         rp1.x=rp2.x=tmp;
    1276.                         rp1.y=sqrt(r*r-tmp*tmp);
    1277.                         rp2.y=-rp1.y;
    1278.                 }
    1279.         }
    1280.         else if(a!=0&&b!=0)
    1281.         {
    1282.                 double delta;
    1283.                 delta=b*b*c*c-(a*a+b*b)*(c*c-a*a*r*r);
    1284.                 if(delta<0)
    1285.                         res=0;
    1286.                 else if(delta==0)
    1287.                 {
    1288.                         res=1;
    1289.                         rp1.y=-b*c/(a*a+b*b);
    1290.                         rp1.x=(-c-b*rp1.y)/a;
    1291.                 }
    1292.                 else
    1293.                 {
    1294.                         res=2;
    1295.                         rp1.y=(-b*c+sqrt(delta))/(a*a+b*b);
    1296.                         rp2.y=(-b*c-sqrt(delta))/(a*a+b*b);
    1297.                         rp1.x=(-c-b*rp1.y)/a;
    1298.                         rp2.x=(-c-b*rp2.y)/a;
    1299.                 }
    1300.         }
    1301.         rp1.x+=p.x;
    1302.         rp1.y+=p.y;
    1303.         rp2.x+=p.x;
    1304.         rp2.y+=p.y;
    1305.         return res;
    1306. }
    1307. //内切圆:
    1308. void incircle(POINT p1,POINT p2,POINT p3,POINT &rp,double &r)
    1309. {
    1310.         double dx31,dy31,dx21,dy21,d31,d21,a1,b1,c1;
    1311.         dx31=p3.x-p1.x;
    1312.         dy31=p3.y-p1.y;
    1313.         dx21=p2.x-p1.x;
    1314.         dy21=p2.y-p1.y;

    1315.         d31=sqrt(dx31*dx31+dy31*dy31);
    1316.         d21=sqrt(dx21*dx21+dy21*dy21);
    1317.         a1=dx31*d21-dx21*d31;
    1318.         b1=dy31*d21-dy21*d31;
    1319.         c1=a1*p1.x+b1*p1.y;

    1320.         double dx32,dy32,dx12,dy12,d32,d12,a2,b2,c2;
    1321.         dx32=p3.x-p2.x;
    1322.         dy32=p3.y-p2.y;
    1323.         dx12=-dx21;
    1324.         dy12=-dy21;

    1325.         d32=sqrt(dx32*dx32+dy32*dy32);
    1326.         d12=d21;
    1327.         a2=dx12*d32-dx32*d12;
    1328.         b2=dy12*d32-dy32*d12;
    1329.         c2=a2*p2.x+b2*p2.y;

    1330.         rp.x=(c1*b2-c2*b1)/(a1*b2-a2*b1);
    1331.         rp.y=(c2*a1-c1*a2)/(a1*b2-a2*b1);
    1332.         r=fabs(dy21*rp.x-dx21*rp.y+dx21*p1.y-dy21*p1.x)/d21;
    1333. }
    1334. //求切点:
    1335. // p---圆心坐标, r---圆半径, sp---圆外一点, rp1,rp2---切点坐标
    1336. void cutpoint(POINT p,double r,POINT sp,POINT &rp1,POINT &rp2)
    1337. {
    1338.         POINT p2;
    1339.         p2.x=(p.x+sp.x)/2;
    1340.         p2.y=(p.y+sp.y)/2;

    1341.         double dx2,dy2,r2;
    1342.         dx2=p2.x-p.x;
    1343.         dy2=p2.y-p.y;
    1344.         r2=sqrt(dx2*dx2+dy2*dy2);
    1345.         c2point(p,r,p2,r2,rp1,rp2);
    1346. }
    1347. //线段的左右旋:
    1348. /* l2在l1的左/右方向(l1为基准线)
    1349. 返回        0        : 重合;
    1350. 返回        1        : 右旋;
    1351. 返回        –1 : 左旋;
    1352. */
    1353. int rotat(LINESEG l1,LINESEG l2)
    1354. {
    1355.         double dx1,dx2,dy1,dy2;
    1356.         dx1=l1.s.x-l1.e.x;
    1357.         dy1=l1.s.y-l1.e.y;
    1358.         dx2=l2.s.x-l2.e.x;
    1359.         dy2=l2.s.y-l2.e.y;

    1360.         double d;
    1361.         d=dx1*dy2-dx2*dy1;
    1362.         if(d==0)
    1363.                 return 0;
    1364.         else if(d>0)
    1365.                 return -1;
    1366.         else
    1367.                 return 1;
    1368. }
    1369. /*
    1370. 公式:
    1371. 球坐标公式:
    1372. 直角坐标为 P(x, y, z) 时,对应的球坐标是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP与Z轴的夹角,范围[0,π];是OP在XOY面上的投影到X轴的旋角,范围[0,2π]  
    1373. 直线的一般方程转化成向量方程:
    1374. ax+by+c=0
    1375. x-x0     y-y0
    1376.    ------ = ------- // (x0,y0)为直线上一点,m,n为向量
    1377. m        n
    1378. 转换关系:
    1379. a=n;b=-m;c=m·y0-n·x0;
    1380. m=-b; n=a;
    1381. 三点平面方程:
    1382. 三点为P1,P2,P3
    1383. 设向量  M1=P2-P1; M2=P3-P1;
    1384. 平面法向量:  M=M1 x M2 ()
    1385. 平面方程:    M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0
    1386. */
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-21 16:29 | 显示全部楼层
    1. #include<cstdlib>
    2. #include<cmath>
    3. #include<cstdio>
    4. #include<algorithm>
    5. #define max(a,b) (((a)>(b))?(a)b))
    6. #define min(a,b) (((a)>(b))?(b)a))
    7. #define sign(x) ((x)>eps?1(x)<-eps?(-1)0)))
    8. using namespace std;
    9. const int MAXN=1000;
    10. const double eps=1e-8,inf=1e50 ,Pi=acos(-1.0);
    11. struct point {
    12.     double x,y;
    13.     point() {}
    14.     point(double _x,double _y) {
    15.         x=_x;
    16.         y=_y;
    17.     }
    18.     point operator-(const point &ne)const {
    19.         return point(x-ne.x,y-ne.y);
    20.     }
    21.     point operator+(const point ne)const {
    22.         return point(x+ne.x,y+ne.y);
    23.     }
    24.     point operator*(const double t)const{
    25.         return point(x*t,y*t);
    26.     }
    27.     point operator/(const double t)const{
    28.         if(sign(t)==0)exit(1);
    29.         return point(x/t,y/t);
    30.     }
    31. };
    32. struct line{
    33.     point a,b;
    34.     line(){}
    35.     line(point _a,point _b){a=_a;b=_b;}
    36. };
    37. struct line2 {
    38.     double a,b,c;
    39.     line2() {}
    40.     line2(double _a,double _b,double _c) {
    41.         a=_a;
    42.         b=_b;
    43.         c=_c;
    44.     }
    45. };
    46. struct circle {
    47.     point o;
    48.     double r;
    49.     circle(){}
    50.     circle(point _o,double _r){
    51.         o=_o;
    52.         r=_r;
    53.     }
    54. };
    55. struct rectangle{
    56.     point a,b,c,d;
    57.     rectangle(){}
    58.     rectangle(point _a,point _b,point _c,point _d){
    59.         a=_a;
    60.         b=_b;
    61.         c=_c;
    62.         d=_d;
    63.     }
    64. };
    65. struct polygon{
    66.     point p[MAXN];
    67.     int n;
    68. };
    69. inline double xmult(point a,point b){
    70.     return a.x*b.y-a.y*b.x;
    71. }
    72. inline double xmult(point o,point a,point b) {
    73.     return (a.x-o.x)*(b.y-o.y)-(b.x-o.x)*(a.y-o.y);
    74. }
    75. inline double xmult(double x1,double y1,double x2,double y2) {
    76.     return x1*y2-x2*y1;
    77. }
    78. inline double dmult(point o,point a,point b) {
    79.     return (a.x-o.x)*(b.x-o.x)+(a.y-o.y)*(b.y-o.y);
    80. }
    81. inline double dmult(point a,point b) {
    82.     return a.x*b.x+a.y*b.y;
    83. }
    84. inline double lenth(point a){
    85.     return sqrt(dmult(a,a));
    86. }
    87. inline double dist(point a,point b){
    88.     return lenth(b-a);
    89. }
    90. inline double dist2(point a,point b){
    91.     return dmult(b-a,b-a);
    92. }
    93. //直线一般式转两点式
    94. line toline(double a,double b,double c) {
    95.     if(sign(b)==0)exit(1);
    96.     point A(0,-c/b),B(-c/a,0);
    97.     return line(A,B);
    98. }
    99. //直线两点式转一般式
    100. line2 toline2(point a,point b) {
    101.     double A=b.y-a.y,B=a.x-b.x,C=-B*a.y-A*a.x;
    102.     return line2(A,B,C);
    103. }
    104. //点p绕o逆时针旋转alpha
    105. point rotate(point o,point p,double alpha) {
    106.     point tp;
    107.     p.x-=o.x;
    108.     p.y-=o.y;
    109.     tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;
    110.     tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;
    111.     return tp;
    112. }
    113. //向量u的倾斜角
    114. double angle(point u) {
    115.     return atan2(u.y,u.x);
    116. }
    117. //oe与os的夹角,夹角正负满足叉积
    118. double angle(point o,point s,point e) {
    119.     point os=s-o,oe=e-o;
    120.     double bot=sqrt(dmult(os,os)*dmult(oe,oe));
    121.     double top=dmult(os,oe);
    122.     double cosfi=top/bot;
    123.     if (cosfi >=  1.0 ) return 0;
    124.     if (cosfi <= -1.0 ) return -Pi;
    125.     double fi=acos(cosfi);
    126.     if(xmult(o,s,e)>0)return fi;
    127.     else return -fi;
    128. }
    129. //p在l上的投影与l关系
    130. double relation(point p,line l) {
    131.     line tl(l.a,p);
    132.     return dmult(tl.b-l.a,l.b-l.a)/dist2(l.a,l.b);
    133. }
    134. //p在l上的垂足
    135. point perpendicular(point p,line l) {
    136.     double r=relation(p,l);
    137.     return l.a+((l.b-l.a)*r);
    138. }
    139. //求点p到线段l的最短距离,并返回线段上距该点最近的点np
    140. double dist_p_to_seg(point p,line l,point &np) {
    141.     double r=relation(p,l);
    142.     if(r<0) {
    143.     np=l.a;
    144.         return dist(p,l.a);
    145.     }
    146.     if(r>1) {
    147.         np=l.b;
    148.         return dist(p,l.b);
    149.     }
    150.     np=perpendicular(p,l);
    151.     return dist(p,np);
    152. }
    153. //求点p到直线l的最短距离
    154. double dist_p_to_line(point p,line l) {
    155.     return abs(xmult(p-l.a,l.b-l.a))/dist(l.a,l.b);
    156. }
    157. //线段之间最短距离
    158. inline double dist_seg_to_seg(line p,line q) {
    159.     return min(min(dist_p_to_seg(p.a,q),dist_p_to_seg(p.b,q)),min(dist_p_to_seg(q.a,p),dist_p_to_seg(q.b,p)));
    160. }
    161. //求矢量线段夹角的余弦
    162. double cosine(line u,line v) {
    163.     point pu=u.b-u.a,pv=v.b-v.a;
    164.     return dmult(pu,pv)/sqrt(dmult(pu,pu)*dmult(pv,pv));
    165. }
    166. //求矢量的夹角的余弦
    167. double cosine(point a,point b){
    168.     return dmult(a,b)/sqrt(dmult(a,a)*dmult(b,b));
    169. }
    170. //求线段夹角
    171. double lsangle(line u,line v) {
    172.     point o(0,0),a=u.b-u.a,b=v.b-v.a;
    173.     return angle(o,a,b);
    174. }
    175. //求直线斜率
    176. double slope(line2 l) {
    177.     if(abs(l.a) < 1e-20)
    178.         return 0;
    179.     if(abs(l.b) < 1e-20)
    180.         return inf;
    181.     return -(l.a/l.b);
    182. }
    183. //直线倾斜角[0,Pi]
    184. double alpha(line2 l) {
    185.     if(abs(l.a)< eps)
    186.         return 0;
    187.     if(abs(l.b)< eps)
    188.         return Pi/2;
    189.     double k=slope(l);
    190.     if(k>0)
    191.         return atan(k);
    192.     else
    193.         return Pi+atan(k);
    194. }
    195. //点关于直线的对称点
    196. point symmetry(line2 l,point p){
    197.    point tp;
    198.    tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
    199.    tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
    200.    return tp;
    201. }
    202. //判多边形是否逆时针
    203. bool is_unclock(polygon pg) {
    204.     int n=pg.n;
    205.     pg.p[n]=pg.p[0];
    206.     double area=0;
    207.     for(int i=0; i<n; i++)
    208.         area+=xmult(pg.p[i].x,pg.p[i].y,pg.p[i+1].x,pg.p[i+1].y);
    209.     return area>-eps;
    210. }
    211. //改变多边形时针顺序
    212. void to_unclock(polygon &pg) {
    213.     for(int i=0,j=pg.n-1; i<j; i++,j--)
    214.         swap(pg.p[i],pg.p[j]);
    215. }

    216. //判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线
    217. int is_convex(point p[],int n) {
    218.     int i,s[3]= {1,1,1};
    219.     for (i=0; i<n&&s[1]|s[2]; i++)
    220.         s[(sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))+3)%3]=0;
    221.     return s[1]|s[2];
    222. }
    223. //判定凸多边形,顶点按顺时针或逆时针给出,不允许相邻边共线
    224. int is_convex_v2(point p[],int n) {
    225.     int i,s[3]= {1,1,1};
    226.     for (i=0; i<n&&s[0]&&s[1]|s[2]; i++)
    227.         s[(sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))+3)%3]=0;
    228.     return s[0]&&s[1]|s[2];
    229. }
    230. //判点在凸多边形内或多边形边上,顶点按顺时针或逆时针给出
    231. int inside_convex(point q,point p[],int n) {
    232.     int i,s[3]= {1,1,1};
    233.     for (i=0; i<n&&s[1]|s[2]; i++)
    234.         s[(sign(xmult(p[(i+1)%n],q,p[i]))+3)%3]=0;
    235.     return s[1]|s[2];
    236. }
    237. //判点在凸多边形内,顶点按顺时针或逆时针给出,在多边形边上返回0
    238. int inside_convex2(point q,point p[],int n) {
    239.     int i,s[3]= {1,1,1};
    240.     for (i=0; i<n&&s[0]&&s[1]|s[2]; i++)
    241.         s[(sign(xmult(p[(i+1)%n],q,p[i]))+3)%3]=0;
    242.     return s[0]&&s[1]|s[2];
    243. }
    244. //判点在线段上
    245. inline int p_on_seg(point a,point p1,point p2) {
    246.     if(fabs(xmult(a,p1,p2))<=eps&&(a.x-p1.x)*(a.x-p2.x)<eps&&(a.y-p1.y)*(a.y-p2.y)<eps)
    247.         return 1;
    248.     return 0;
    249. }
    250. //判点在线段端点左方
    251. inline int p_on_segvex(point s,point p) {
    252.     return fabs(p.y-s.y)<eps&&(p.x<=s.x+eps);
    253. }
    254. //判线段相交 <=:不规范相交
    255. inline int seg_inter(line s,line p) {
    256.     double minx1=min(s.a.x,s.b.x),maxx1=max(s.a.x,s.b.x);
    257.     double minx2=min(p.a.x,p.b.x),maxx2=max(p.a.x,p.b.x);
    258.     double miny1=min(s.a.y,s.b.y),maxy1=max(s.a.y,s.b.y);
    259.     double miny2=min(p.a.y,p.b.y),maxy2=max(p.a.y,p.b.y);
    260.     if((minx1>maxx2+eps)||(minx2>maxx1+eps)||
    261.             (miny1>maxy2+eps)||(miny2>maxy1+eps))
    262.         return 0;
    263.     else
    264.         return sign(xmult(s.a,s.b,p.a)*xmult(s.a,s.b,p.b))<=0&&
    265.                sign(xmult(p.a,p.b,s.a)*xmult(p.a,p.b,s.b))<=0;
    266. }
    267. //判点在多边形内部
    268. inline int p_in_polygon(point a,point p[],int n) {
    269.     int count = 0;
    270.     line s,ps;
    271.     ps.a = a,ps.b = a;
    272.     ps.b.x = inf;
    273.     for(int i = 0; i < n; i++) {
    274.         s.a = p[i];
    275.         if(i + 1 < n)s.b = p[i+1];
    276.         else s.b = p[0];
    277.         if (s.a.y > s.b.y)swap(s.a,s.b);
    278.         if (p_on_seg(a,s.a,s.b))return 2;
    279.         if ((fabs(s.a.y-s.b.y)>eps)) {
    280.             if (p_on_segvex(s.b,a)) count++;
    281.             else if (seg_inter(ps,s))count++;
    282.         }
    283.     }
    284.     if (count%2)return 1;
    285.     return 0;
    286. }
    287. //多边形内部最长线段
    288. point stk[MAXN];
    289. double seg_max_len(line u,polygon &pg){
    290.     double ans=0.0,tmp=inf;
    291.     pg.p[pg.n]=pg.p[0];
    292.     int n=pg.n,top=0;
    293.     for(int i=0;i<n;i++){
    294.         line v(pg.p[i],pg.p[i+1]);
    295.         int s1=sign(xmult(u.a,u.b,v.a)),s2=sign(xmult(u.a,u.b,v.b));
    296.         if(s1*s2<=0&&(s1!=0||s2!=0)){
    297.             stk[top++]=line_intersection(u,v);
    298.         }
    299.     }
    300.     stk[top++]=u.a;
    301.     stk[top++]=u.b;
    302.     sort(stk,stk+top);
    303.     top=unique(stk,stk+top)-stk;
    304.     point mp,lp=stk[0];
    305.     for(int i=1;i<top;i++){
    306.         mp=(lp+stk[i])*0.5;
    307.         if(!p_in_polygon(mp,pg))lp=stk[i];
    308.         ans=max(ans,dist2(lp,stk[i]));
    309.     }
    310.     return sqrt(ans);
    311. }
    312. double maxlenth(polygon &pg){
    313.     double ans=0.0;
    314.     for(int i=0;i<pg.n-1;i++){
    315.         for(int j=i+1;j<pg.n;j++){
    316.             ans=max(ans,seg_max_len(line(pg.p[i],pg.p[j]),pg));
    317.             ans=max(ans,seg_max_len(line(pg.p[j],pg.p[i]),pg));
    318.         }
    319.     }
    320.     return ans;
    321. }
    322. //凸包对踵点长度
    323. double opposite_lenth(polygon &pg){
    324.     double ans=inf;
    325.     int a,b,c;
    326.     pg.p[pg.n]=pg.p[0];
    327.     for(a=0,b=1,c=2;a<pg.n;a++,b++){
    328.         while(getarea(pg.p[a],pg.p[b],pg.p[c])<getarea(pg.p[a],pg.p[b],pg.p[(c+1)%pg.n]))c=(c+1)%pg.n;
    329.         ans=min(ans,dist_p_to_line(pg.p[c],line(pg.p[a],pg.p[b])));
    330.     }
    331.     return ans;
    332. }
    333. //判p1,p2是否在l1,l2两侧
    334. inline int opposite_side(point p1,point p2,point l1,point l2) {
    335.     return xmult(l1,l2,p1)*xmult(l1,l2,p2)<-eps;
    336. }
    337. //判线段在任意多边形内,顶点按顺时针或逆时针给出,与边界相交返回1
    338. int seg_in_polygon(point l1,point l2,point p[],int n) {
    339.     point t[MAXN],tt;
    340.     int i,j,k=0;
    341.     if (!p_in_polygon(l1,p,n)||!p_in_polygon(l2,p,n))//不在内部
    342.         return 0;
    343.     for (i=0; i<n; i++)
    344.         if (opposite_side(l1,l2,p[i],p[(i+1)%n])&&opposite_side(p[i],p[(i+1)%n],l1,l2))
    345.             return 0;
    346.         else if (p_on_seg(l1,p[i],p[(i+1)%n]))
    347.             t[k++]=l1;
    348.         else if (p_on_seg(l2,p[i],p[(i+1)%n]))
    349.             t[k++]=l2;
    350.         else if (p_on_seg(p[i],l1,l2))
    351.             t[k++]=p[i];
    352.     for (i=0; i<k; i++)
    353.         for (j=i+1; j<k; j++) {
    354.             tt.x=(t[i].x+t[j].x)/2;
    355.             tt.y=(t[i].y+t[j].y)/2;
    356.             if (!p_in_polygon(tt,p,n))
    357.                 return 0;
    358.         }
    359.     return 1;
    360. }
    361. //求直线交点,必须存在交点,或者预判断【解析几何方法】
    362. point line_intersection(line u,line v) {
    363.     double a1=u.b.y-u.a.y,b1=u.a.x-u.b.x;
    364.     double c1=u.b.y*(-b1)-u.b.x*a1;
    365.     double a2=v.b.y-v.a.y,b2=v.a.x-v.b.x;
    366.     double c2=v.b.y*(-b2)-v.b.x*a2;
    367.     double D=xmult(a1,b1,a2,b2);
    368.     return point(xmult(b1,c1,b2,c2)/D,xmult(c1,a1,c2,a2)/D);
    369. }
    370. //求线段交点,必须存在交点,或者预判断【平面几何方法】
    371. point line_intersection2(line u,line v) {
    372.     point ret=u.a;
    373.     double t=xmult(u.a-v.a,v.b-v.a)/xmult(u.b-u.a,v.b-v.a);
    374.     t=fabs(t);
    375.     ret.x+=(u.b.x-u.a.x)*t;
    376.     ret.y+=(u.b.y-u.a.y)*t;
    377.     return ret;
    378. }
    379. //三角形重心
    380. point barycenter(point a,point b,point c){
    381.     return (a+b+c)/3.0;
    382. }
    383. //多边形重心
    384. point barycenter(point p[],int n) {
    385.     point ret,t;
    386.     double t1=0,t2;
    387.     int i;
    388.     ret.x=ret.y=0;
    389.     for (i=1; i<n-1; i++)
    390.         if (fabs(t2=xmult(p[i+1],p[0],p[i]))>eps) {
    391.             t=barycenter(p[0],p[i],p[i+1]);
    392.             ret.x+=t.x*t2;
    393.             ret.y+=t.y*t2;
    394.             t1+=t2;
    395.         }
    396.     if (fabs(t1)>eps)
    397.         ret.x/=t1,ret.y/=t1;
    398.     return ret;
    399. }
    400. //求多边形面积
    401. inline double getarea(point pg[],int n) {
    402.     double area=0;
    403.     pg[n]=pg[0];
    404.     for(int i=0; i<n; i++)
    405.         area+=xmult(pg[i].x,pg[i].y,pg[i+1].x,pg[i+1].y);
    406.     return fabs(area)/2.0;
    407. }
    408. //解方程ax^2+bx+c=0
    409. int equaltion(double a,double b,double c,double &x1,double &x2) {
    410.     double der=b*b-4*a*c;
    411.     switch(sign(der)) {
    412.     case -1:
    413.         return 0;
    414.     case 0:
    415.         x1=-b/(2*a);
    416.         return 1;
    417.     case 1:
    418.         der=sqrt(der);
    419.         x1=(-b-der)/(2*a);
    420.         x2=(-b+der)/(2*a);
    421.         return 2;
    422.     }
    423. }
    424. //线段与圆交点
    425. int line_circle_intersection(line u,circle c,point &p1,point &p2) {
    426.     double dis=lenth(u.b-u.a);
    427.     point d=(u.b-u.a)/dis;
    428.     point E=c.o-u.a;
    429.     double a=dmult(E,d);
    430.     double a2=a*a;
    431.     double e2=dmult(E,E);
    432.     double r2=c.r*c.r;
    433.     if((r2-e2+a2)<0){
    434.         return 0;
    435.     } else {
    436.         double f=sqrt(r2 - e2 + a2);
    437.         double t=a-f;
    438.         int cnt=0;
    439.         if(t>-eps&&t-dis<eps){//去掉后面变成射线
    440.             p1=u.a+(d*t);
    441.             cnt++;
    442.         }
    443.         t=a+f;
    444.         if(t>-eps&&t-dis<eps) {
    445.             p2=u.a+(d*t);
    446.             cnt++;
    447.         }
    448.         return cnt;
    449.     }
    450. }
    451. //给出在任意多边形内部的一个点
    452. point a_point_in_polygon(polygon pg) {
    453.     point v,a,b,r;
    454.     int i,index;
    455.     v=pg.p[0];
    456.     index=0;
    457.     for(i=1; i<pg.n; i++) {
    458.         if(pg.p[i].y<v.y) {
    459.             v=pg.p[i];
    460.             index=i;
    461.         }
    462.     }
    463.     a=pg.p[(index-1+pg.n)%pg.n];
    464.     b=pg.p[(index+1)%pg.n];
    465.     point q;
    466.     polygon tri;
    467.     tri.n=3;
    468.     tri.p[0]=a;
    469.     tri.p[1]=v;
    470.     tri.p[2]=b;
    471.     double md=inf;
    472.     int in1=index;
    473.     bool bin=false;
    474.     for(i=0; i<pg.n; i++) {
    475.         if(i == index)continue;
    476.         if(i == (index-1+pg.n)%pg.n)continue;
    477.         if(i == (index+1)%pg.n)continue;
    478.         if(!inside_convex2(pg.p[i],tri.p,3))continue;

    479.         bin=true;
    480.         if(dist(v,pg.p[i])<md) {
    481.             q=pg.p[i];
    482.             md=dist(v,q);
    483.         }
    484.     }
    485.     if(!bin) {
    486.         r.x=(a.x+b.x)/2;
    487.         r.y=(a.y+b.y)/2;
    488.         return r;
    489.     }
    490.     r.x=(v.x+q.x)/2;
    491.     r.y=(v.y+q.y)/2;
    492.     return r;
    493. }
    494. //求在多边形外面的点到凸包的切点
    495. void p_cut_polygon(point p,polygon pg,point &rp,point &lp) {
    496.     line ep,en;
    497.     bool blp,bln;
    498.     rp=pg.p[0];
    499.     lp=pg.p[0];
    500.     for(int i=1; i<pg.n; i++) {
    501.         ep.a=pg.p[(i+pg.n-1)%pg.n];
    502.         ep.b=pg.p[i];
    503.         en.a=pg.p[i];
    504.         en.b=pg.p[(i+1)%pg.n];
    505.         blp=xmult(ep.b-ep.a,p-ep.a)>=0;
    506.         bln=xmult(en.b-en.a,p-en.a)>=0;
    507.         if(!blp&&bln) {
    508.             if(xmult(pg.p[i]-p,rp-p)>0)
    509.                 rp=pg.p[i];
    510.         }
    511.         if(blp&&!bln) {
    512.             if(xmult(lp-p,pg.p[i]-p)>0)
    513.                 lp=pg.p[i];
    514.         }
    515.     }
    516.     return ;
    517. }
    518. //判断点p在圆c内
    519. bool p_in_circle(point p,circle c) {
    520.     return c.r*c.r>dist2(p,c.o);
    521. }
    522. //求矩形第4个点
    523. point rect4th(point a,point b,point c) {
    524.     point d;
    525.     if(abs(dmult(a-c,b-c))<eps) {
    526.         d=a+b-c;
    527.     }
    528.     if(abs(dmult(a-b,c-b))<eps) {
    529.         d=a+c-b;
    530.     }
    531.     if(abs(dmult(c-a,b-a))<eps) {
    532.         d=c+b-a;
    533.     }
    534.     return d;
    535. }
    536. //判两圆关系
    537. int CircleRelation(circle c1,circle c2){
    538.     double d=lenth(c1.o-c2.o);
    539.     if( fabs(d-c1.r-c2.r) < eps )
    540.         return 2;//外切
    541.     if( fabs(d-fabs(c1.r-c2.r)) < eps )
    542.         return 4;//内切
    543.     if( d > c1.r+c2.r )
    544.         return 1;//相离
    545.     if( d < fabs(c1.r-c2.r) )
    546.         return 5;//内含
    547.     if( fabs(c1.r-c2.r) < d && d < c1.r+c2.r)
    548.         return 3;//相交
    549.     return 0; //error!
    550. }
    551. //判圆与矩形关系,矩形水平
    552. bool Circle_In_Rec(circle c,rectangle r) {
    553.     if( r.a.x < c.o.x && c.o.x < r.b.x && r.c.y < c.o.y && c.o.y < r.b.y ) {
    554.         line line1(r.a, r.b);
    555.         line line2(r.b, r.c);
    556.         line line3(r.c, r.d);
    557.         line line4(r.d, r.a);
    558.         if(c.r<dist_p_to_line(c.o,line1)&&c.r<dist_p_to_line(c.o,line2)&&c.r<dist_p_to_line(c.o,line3)&&c.r<dist_p_to_line(c.o,line4))
    559.             return true;
    560.     }
    561.     return false;
    562. }
    563. //射线关于平面的反射
    564. void reflect(line2 u,line2 v,line2 &l){
    565.     double n,m;
    566.     double tpb,tpa;
    567.     tpb=u.b*v.b+u.a*v.a;
    568.     tpa=v.a*u.b-u.a*v.b;
    569.     m=(tpb*u.b+tpa*u.a)/(u.b*u.b+u.a*u.a);
    570.     n=(tpa*u.b-tpb*u.a)/(u.b*u.b+u.a*u.a);
    571.     if(fabs(u.a*v.b-v.a*u.b)<1e-20) {
    572.         l.a=v.a;
    573.         l.b=v.b;
    574.         l.c=v.c;
    575.         return;
    576.     }
    577.     double xx,yy; //(xx,yy)是入射线与镜面的交点。
    578.     xx=(u.b*v.c-v.b*u.c)/(u.a*v.b-v.a*u.b);
    579.     yy=(v.a*u.c-u.a*v.c)/(u.a*v.b-v.a*u.b);
    580.     l.a=n;
    581.     l.b=-m;
    582.     l.c=m*yy-xx*n;
    583. }
    584. //两圆交点(预判断不相交情况)
    585. void c2point(circle c1,circle c2,point &rp1,point &rp2) {
    586.     double a,b,r;
    587.     a=c2.o.x-c1.o.x;
    588.     b=c2.o.y-c1.o.y;
    589.     r=(a*a+b*b+c1.r*c1.r-c2.r*c2.r)/2;
    590.     if(a==0&&b!=0) {
    591.         rp1.y=rp2.y=r/b;
    592.         rp1.x=sqrt(c1.r*c1.r-rp1.y*rp1.y);
    593.         rp2.x=-rp1.x;
    594.     } else if(a!=0&&b==0) {
    595.         rp1.x=rp2.x=r/a;
    596.         rp1.y=sqrt(c1.r*c1.r-rp1.x*rp2.x);
    597.         rp2.y=-rp1.y;
    598.     } else if(a!=0&&b!=0) {
    599.         double delta;
    600.         delta=b*b*r*r-(a*a+b*b)*(r*r-c1.r*c1.r*a*a);
    601.         rp1.y=(b*r+sqrt(delta))/(a*a+b*b);
    602.         rp2.y=(b*r-sqrt(delta))/(a*a+b*b);
    603.         rp1.x=(r-b*rp1.y)/a;
    604.         rp2.x=(r-b*rp2.y)/a;
    605.     }
    606.     rp1=rp1+c1.o;
    607.     rp2=rp2+c1.o;
    608. }
    609. //圆外一点引圆的切线
    610. void cutpoint(circle c,point sp,point &rp1,point &rp2) {
    611.     circle c2;
    612.     c2.o=(c.o+sp)/2.0;
    613.     c2.r=lenth(c2.o-sp);
    614.     c2point(c,c2,rp1,rp2);
    615. }
    616. //圆c1上,与c2的外切点
    617. void c2cuto(circle c1, circle c2, point &p1, point &p2) {
    618.     double d = dist(c1.o, c2.o), dr = c1.r - c2.r;
    619.     double b = acos(dr / d);
    620.     double a = angle(c2.o-c1.o);
    621.     double a1 = a - b, a2 = a + b;
    622.     p1=point(cos(a1) * c1.r, sin(a1) * c1.r) + c1.o;
    623.     p2=point(cos(a2) * c1.r, sin(a2) * c1.r) + c1.o;
    624. }
    625. //圆c1上,与c2的内切点
    626. void c2cuti(circle c1,circle c2,point &p1,point &p2){
    627.     point dr=c2.o-c1.o;
    628.     dr=dr/lenth(dr);
    629.     point a=c1.o-(dr*c1.r),b=c1.o+(dr*c1.r);
    630.     point c=c2.o-(dr*c2.r),d=c2.o+(dr*c2.r);
    631.     circle E((a+c)/2.0,lenth(c-a)/2.0),F((b+d)/2.0,lenth(d-b)/2.0);
    632.     point q1,q2;
    633.     c2point(E,F,q1,q2);
    634.     point L=line_intersection2(line(c1.o,c2.o),line(q1,q2));
    635.     circle c3((c1.o+L)/2.0,lenth(L-c1.o)/2.0);
    636.     c2point(c1,c3,p1,p2);
    637. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-21 16:30 | 显示全部楼层

    1. #include"string.h"
    2. #include"stdio.h"
    3. #include"iostream"
    4. #include"algorithm"
    5. #include"queue"
    6. #include"stack"
    7. #define M 100009
    8. #define N 100009
    9. #include"stdlib.h"
    10. #include"math.h"
    11. #define inf 10000000000000000LL
    12. #define INF 0x3f3f3f3f
    13. #define PI acos(-1.0)
    14. #define eps 1e-10
    15. using namespace std;
    16. struct node
    17. {
    18.     double x,y;
    19.     node(){}
    20.     node(double _x,double _y):x(_x),y(_y){}
    21.     node operator +(node p)//向量加法
    22.     {
    23.         return node(x+p.x,y+p.y);
    24.     }
    25.     node operator -(node p)//向量减法
    26.     {
    27.         return node(x-p.x,y-p.y);
    28.     }
    29.     double operator *(node p)//向量叉乘
    30.     {
    31.         return x*p.y-y*p.x;
    32.     }
    33.     double operator ^(node p)//向量点乘
    34.     {
    35.         return x*p.x+y*p.y;
    36.     }
    37.     node operator /(double p)//向量除法
    38.     {
    39.         return node(x/p,y/p);
    40.     }
    41.     node operator *(double p)//向量乘法
    42.     {
    43.         return node(x*p,y*p);
    44.     }
    45. }p[M],q[M];
    46. int cnt,n;
    47. double max(double x,double y)
    48. {
    49.     return x>y?x:y;
    50. }
    51. double min(double x,double y)
    52. {
    53.     return x<y?x:y;
    54. }
    55. double cross(node a,node b,node c)//叉积
    56. {
    57.     return (b-a)*(c-a);
    58. }
    59. double dot(node a,node b,node c)//点积
    60. {
    61.     return (b-a)^(c-a);
    62. }
    63. double len(node a)//向量长吨
    64. {
    65.     return sqrt(a^a);
    66. }
    67. double dis(node a,node b)//两点距离
    68. {
    69.     return len(b-a);
    70. }
    71. int cmp(node a,node b)//极角排序
    72. {
    73.     double temp=cross(p[0],a,b);//逆时针排序
    74.     if(temp>0)
    75.         return 1;
    76.     else if(fabs(temp)<eps&&dis(p[0],a)<dis(p[0],b))//角度相同则按照距离排序
    77.         return 1;
    78.     else
    79.         return 0;
    80. }
    81. void input(int n)//输入
    82. {
    83.     for(int i=0;i<n;i++)
    84.         scanf("%lf%lf",&p[i].x,&p[i].y);
    85. }
    86. void sort_point(int n)//凸包点集的输入即排序
    87. {
    88.     int i,k;
    89.     node start;
    90.     start=p[0];
    91.     k=0;
    92.     for(i=1;i<n;i++)
    93.     {
    94.         if((start.y>p[i].y)||(fabs(start.y-p[i].y)<eps&&start.x>p[i].x))
    95.         {
    96.             start=p[i];
    97.             k=i;
    98.         }
    99.     }
    100.     p[k]=p[0];
    101.     p[0]=start;
    102.     sort(p+1,p+n,cmp);
    103. }
    104. void be_weight(int val)
    105. {
    106.     int temp=val;
    107.     n=1;
    108.     for(int i=1;i<temp;i++)
    109.     {
    110.         if(fabs(p[i-1].x-p[i].x)<eps&&fabs(p[i-1].y-p[i].y)<eps)
    111.             continue;
    112.         p[n++]=p[i];
    113.     }
    114. }
    115. void Convex_hull(int n)//求凸包凸包上的点存在q中
    116. {
    117.     int i;
    118.     if(n==1)
    119.     {
    120.         q[0]=p[0];
    121.         cnt=1;
    122.     }
    123.     else if(n==2)
    124.     {
    125.         q[0]=p[0];
    126.         q[1]=p[1];
    127.         q[2]=p[0];
    128.         cnt=2;
    129.     }
    130.     else if(n>=3)
    131.     {
    132.         q[0]=p[n-1];
    133.         q[1]=p[0];
    134.         q[2]=p[1];
    135.         cnt=2;
    136.         for(i=2;i<n;i++)
    137.         {
    138.             while(cross(q[cnt-1],q[cnt],p[i])<0)
    139.                 cnt--;
    140.             q[++cnt]=p[i];
    141.         }
    142.     }
    143. }
    144. double Perimeter(int cnt)//凸包周长
    145. {
    146.     double sum=0;
    147.     for(int i=1;i<=cnt;i++)
    148.         sum+=dis(q[i-1],q[i]);
    149.     return sum;
    150. }
    151. double Area(int cnt)//凸包面积
    152. {
    153.     double sum=0;
    154.     node p(0,0);
    155.     for(int i=1;i<=cnt;i++)
    156.         sum+=cross(p,q[i-1],q[i]);
    157.     return fabs(sum/2.0);
    158. }
    159. node barycenter_cur(int n)//原多边形的重心
    160. {
    161.     double sum=0;
    162.     node ret(0.0,0.0);
    163.     for(int i=2;i<n;i++)
    164.     {
    165.         double area=cross(p[0],p[i-1],p[i]);
    166.         sum+=area;
    167.         ret=ret+(p[0]+p[i-1]+p[i])/3.0*area;
    168.     }
    169.     ret=ret/sum;
    170.     return ret;
    171. }
    172. node barycenter_now(int cnt)//凸包的重心
    173. {
    174.     double sum=0;
    175.     node ret(0.0,0.0);
    176.     for(int i=2;i<cnt;i++)
    177.     {
    178.         double area=cross(q[0],q[i-1],q[i]);
    179.         sum+=area;
    180.         ret=ret+(q[0]+q[i-1]+q[i])/3.0*area;
    181.     }
    182.     ret=ret/sum;
    183.     return ret;
    184. }
    185. double Diameter(int cnt)//旋转卡壳法求凸包的直径即最大的点对距离
    186. {
    187.     double maxi=0;
    188.     int j=1;
    189.     for(int i=1;i<=cnt;i++)
    190.     {
    191.         while(fabs(cross(q[i-1],q[i],q[(j+1)%cnt]))>fabs(cross(q[i-1],q[i],q[j%cnt])))
    192.             j++;
    193.         maxi=max(maxi,dis(q[i-1],q[j%cnt]));
    194.         maxi=max(maxi,dis(q[i],q[j%cnt]));
    195.     }
    196.     return maxi;
    197. }
    198. double Max_triangle(int cnt)//旋转卡壳法求面积最大的三角形
    199. {
    200.     double maxi=0;
    201.     int j=1;
    202.     int k=1;
    203.     for(int i=1;i<=cnt;i++)
    204.     {
    205.         while(fabs(cross(q[i-1],q[j%cnt],q[(k+1)%cnt]))>fabs(cross(q[i-1],q[j%cnt],q[k%cnt])))
    206.             k++;
    207.         maxi=max(maxi,fabs(cross(q[i-1],q[j%cnt],q[k%cnt])));
    208.         while(fabs(cross(q[i-1],q[(j+1)%cnt],q[k%cnt]))>fabs(cross(q[i-1],q[j%cnt],q[k%cnt])))
    209.             j++;
    210.         maxi=max(maxi,fabs(cross(q[i-1],q[j%cnt],q[k%cnt])));
    211.     }
    212.     return maxi/2.0;
    213.     /*其思路是这样的,定点I,p,q,先I,p固定,让q旋转找到最大的面积三角形,之后,I,q固定,p旋转,
    214.         找到最大的三角形面积,比较记录.然后i++;直到i遍历所有顶点.所求出来的三角形就是面积
    215.         最大.这里的旋转卡壳思想就是固定,旋转.这样的.显然i++后,p,q两点不需要再从i+1,i+2开始,这
    216.         个好形容,对p,q进行取模运算的时候,注意自己的SP栈指针多大.*/
    217. }
    218. void Min_rectangle(int cnt)//旋转卡壳法求面积和周长最小的环绕矩形
    219. {
    220.     if(cnt<=2)//输出时注意的地方*****
    221.     {
    222.         if(cnt==1)
    223.             printf("%.2lf %.2lf\n",0.0,0.0);
    224.         else
    225.             printf("%.2lf %.2lf\n",0.0,2*dis(p[0],p[1]));
    226.         return;
    227.     }
    228.     double S=inf,C=inf;
    229.     int j,k,r;
    230.     double h,w;
    231.     j=k=r=1;
    232.     for(int i=1;i<=cnt;i++)
    233.     {
    234.         double L=dis(q[i-1],q[i]);
    235.         while(fabs(cross(q[i-1],q[i],q[(j+1)%cnt]))>fabs(cross(q[i-1],q[i],q[j%cnt])))
    236.             j++;
    237.         h=fabs(cross(q[i-1],q[i],q[j%cnt]))/L;
    238.         while(dot(q[i-1],q[i],q[(k+1)%cnt])>dot(q[i-1],q[i],q[k%cnt]))
    239.             k++;
    240.         if(i==1)
    241.             r=k;
    242.         while(dot(q[i-1],q[i],q[(r+1)%cnt])<=dot(q[i-1],q[i],q[r%cnt]))
    243.             r++;
    244.         w=(dot(q[i-1],q[i],q[k%cnt])-dot(q[i-1],q[i],q[r%cnt]))/L;
    245.         S=min(S,w*h);
    246.         C=min(C,(w+h)*2);
    247.     }
    248.     printf("%.2lf ",S);//输出时注意的地方*****
    249.     printf("%.2lf\n",C);
    250. }
    251. int main()
    252. {
    253.     while(scanf("%d",&n),n)
    254.     {
    255.         input(n);//输入
    256.         sort_point(n);//极角排序
    257.         be_weight(n);//去重
    258.         Convex_hull(n);//求凸包
    259.         Min_rectangle(cnt);
    260.     }
    261.     return 0;
    262. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-21 16:40 | 显示全部楼层
    1. 判断3D点,3D向量,2D点,2D向量,实数是否相等的ARX函数代码。 [复制链接]
    2. admin

    3. TA的每日心情
    4.         开心
    5. 2021-1-21 00:14
    6. 签到天数: 1397 天

    7. 连续签到: 16 天

    8. [LV.10]以坛为家III

    9. 2854

    10. 主题       
    11. 3730

    12. 帖子       
    13. 214748万

    14. 积分
    15. 管理员

    16. Rank: 9Rank: 9Rank: 9

    17. 积分2147483647
    18. IP编辑禁止帖子清理       
    19. 电梯直达
    20. 跳转到指定楼层楼主
    21. 发表于 2020-5-20 15:52 | 只看该作者 回帖奖励
    22. [pcode=cpp,true]
    23. static BOOL                fuzzyEqual(const AcGePoint3d& pt1, const AcGePoint3d& pt2);
    24. static BOOL                fuzzyEqual(const AcGeVector3d& vec1, const AcGeVector3d& vec2);
    25. static BOOL                fuzzyEqual(const AcGePoint2d& pt1, const AcGePoint2d& pt2);
    26. static BOOL                fuzzyEqual(const AcGeVector2d& vec1, const AcGeVector2d& vec2);
    27. static BOOL                fuzzyEqual(double x1, double x2);
    28. [/pcode]

    29. [pcode=cpp,true]
    30. /****************************************************************************
    31. **
    32. **        XdGeUtils::fuzzyEqual (3D POINTS)
    33. **
    34. **        **jma
    35. **
    36. *************************************/

    37. BOOL
    38. XdGeUtils::fuzzyEqual(const AcGePoint3d& pt1, const AcGePoint3d& pt2)
    39. {
    40.         return((pt1 - pt2).length() <= AcGeContext::gTol.equalPoint());
    41. }

    42. /****************************************************************************
    43. **
    44. **        XdGeUtils::fuzzyEqual (3D VECTORS)
    45. **
    46. **        **jma
    47. **
    48. *************************************/

    49. BOOL
    50. XdGeUtils::fuzzyEqual(const AcGeVector3d& vec1, const AcGeVector3d& vec2)
    51. {
    52.         return((vec1 - vec2).length() <= AcGeContext::gTol.equalVector());
    53. }

    54. /****************************************************************************
    55. **
    56. **        XdGeUtils::fuzzyEqual (2D POINTS)
    57. **
    58. **        **jma
    59. **
    60. *************************************/

    61. BOOL
    62. XdGeUtils::fuzzyEqual(const AcGePoint2d& pt1, const AcGePoint2d& pt2)
    63. {
    64.         return((pt1 - pt2).length() <= AcGeContext::gTol.equalPoint());
    65. }

    66. /****************************************************************************
    67. **
    68. **        XdGeUtils::fuzzyEqual (2D VECTORS)
    69. **
    70. **        **jma
    71. **
    72. *************************************/

    73. BOOL
    74. XdGeUtils::fuzzyEqual(const AcGeVector2d& vec1, const AcGeVector2d& vec2)
    75. {
    76.         return((vec1 - vec2).length() <= AcGeContext::gTol.equalVector());
    77. }

    78. /****************************************************************************
    79. **
    80. **        XdGeUtils::fuzzyEqual (DOUBLES)
    81. **
    82. **        **jma
    83. **
    84. *************************************/

    85. BOOL
    86. XdGeUtils::fuzzyEqual(double x1, double x2)
    87. {
    88.         return(fabs(x1 - x2) <= AcGeContext::gTol.equalPoint());
    89. }
    90. [/pcode]
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-21 16:42 | 显示全部楼层
    1. /点表比较函数,配合std::sort使用
    2. //点按xyz排序,小写字母表示从小到大排序,大写则从大到小排序,在前的字母先排序
    3. int cmp_xyz (AcGePoint3d p1,AcGePoint3d p2)
    4. {
    5.         if (p1.x != p2.x)
    6.                 return p1.x < p2.x?1:0;
    7.         else if (p1.y != p2.y)
    8.                 return p1.y < p2.y?1:0;
    9.         else
    10.                 return p1.z < p2.z?1:0;
    11. }
    12. int cmp_xYz (AcGePoint3d p1,AcGePoint3d p2)
    13. {
    14.         if (p1.x != p2.x)
    15.                 return p1.x < p2.x?1:0;
    16.         else if (p1.y != p2.y)
    17.                 return p1.y > p2.y?1:0;
    18.         else
    19.                 return p1.z < p2.z?1:0;
    20. }
    21. int cmp_xYZ (AcGePoint3d p1,AcGePoint3d p2)
    22. {
    23.         if (p1.x != p2.x)
    24.                 return p1.x < p2.x?1:0;
    25.         else if (p1.y != p2.y)
    26.                 return p1.y > p2.y?1:0;
    27.         else
    28.                 return p1.z > p2.z?1:0;
    29. }
    30. int cmp_xyZ (AcGePoint3d p1,AcGePoint3d p2)
    31. {
    32.         if (p1.x != p2.x)
    33.                 return p1.x < p2.x?1:0;
    34.         else if (p1.y != p2.y)
    35.                 return p1.y < p2.y?1:0;
    36.         else
    37.                 return p1.z > p2.z?1:0;
    38. }
    39. int cmp_XYZ (AcGePoint3d p1,AcGePoint3d p2)
    40. {
    41.         if (p1.x != p2.x)
    42.                 return p1.x > p2.x?1:0;
    43.         else if (p1.y != p2.y)
    44.                 return p1.y > p2.y?1:0;
    45.         else
    46.                 return p1.z > p2.z?1:0;
    47. }
    48. int cmp_XyZ (AcGePoint3d p1,AcGePoint3d p2)
    49. {
    50.         if (p1.x != p2.x)
    51.                 return p1.x > p2.x?1:0;
    52.         else if (p1.y != p2.y)
    53.                 return p1.y < p2.y?1:0;
    54.         else
    55.                 return p1.z > p2.z?1:0;
    56. }

    57. int cmp_Xyz (AcGePoint3d p1,AcGePoint3d p2)
    58. {
    59.         if (p1.x != p2.x)
    60.                 return p1.x > p2.x?1:0;
    61.         else if (p1.y != p2.y)
    62.                 return p1.y < p2.y?1:0;
    63.         else
    64.                 return p1.z < p2.z?1:0;
    65. }

    66. int cmp_XYz (AcGePoint3d p1,AcGePoint3d p2)
    67. {
    68.         if (p1.x != p2.x)
    69.                 return p1.x > p2.x?1:0;
    70.         else if (p1.y != p2.y)
    71.                 return p1.y > p2.y?1:0;
    72.         else
    73.                 return p1.z < p2.z?1:0;
    74. }

    75. int cmp_xy (AcGePoint3d p1,AcGePoint3d p2)
    76. {
    77.         if (p1.x != p2.x)
    78.                 return p1.x < p2.x?1:0;
    79.         else
    80.                 return p1.y < p2.y?1:0;        
    81. }
    82. int cmp_xY (AcGePoint3d p1,AcGePoint3d p2)
    83. {
    84.         if (p1.x != p2.x)
    85.                 return p1.x < p2.x?1:0;
    86.         else
    87.                 return p1.y > p2.y?1:0;        
    88. }
    89. int cmp_Xy (AcGePoint3d p1,AcGePoint3d p2)
    90. {
    91.         if (p1.x != p2.x)
    92.                 return p1.x > p2.x?1:0;
    93.         else
    94.                 return p1.y < p2.y?1:0;        
    95. }
    96. int cmp_XY (AcGePoint3d p1,AcGePoint3d p2)
    97. {
    98.         if (p1.x != p2.x)
    99.                 return p1.x > p2.x?1:0;
    100.         else
    101.                 return p1.y > p2.y?1:0;        
    102. }

    103. int cmp_x(AcGePoint3d p1,AcGePoint3d p2)
    104. {
    105.         return p1.x < p2.x?1:0;
    106. }
    107. int cmp_X(AcGePoint3d p1,AcGePoint3d p2)
    108. {
    109.         return p1.x > p2.x?1:0;
    110. }
    111. int cmp_y(AcGePoint3d p1,AcGePoint3d p2)
    112. {
    113.         return p1.y < p2.y?1:0;
    114. }
    115. int cmp_Y(AcGePoint3d p1,AcGePoint3d p2)
    116. {
    117.         return p1.y > p2.y?1:0;
    118. }
    119. int cmp_z(AcGePoint3d p1,AcGePoint3d p2)
    120. {
    121.         return p1.z < p2.z?1:0;
    122. }
    123. int cmp_Z(AcGePoint3d p1,AcGePoint3d p2)
    124. {
    125.         return p1.z > p2.z?1:0;
    126. }
    127. int cmp_YX (AcGePoint3d p1,AcGePoint3d p2)
    128. {
    129.         if (p1.y != p2.y)
    130.                 return p1.y > p2.y?1:0;
    131.         else
    132.                 return p1.x > p2.x?1:0;        
    133. }
    134. int cmp_Yx (AcGePoint3d p1,AcGePoint3d p2)
    135. {
    136.         if (p1.y != p2.y)
    137.                 return p1.y > p2.y?1:0;
    138.         else
    139.                 return p1.x < p2.x?1:0;        
    140. }
    141. int cmp_yx (AcGePoint3d p1,AcGePoint3d p2)
    142. {
    143.         if (p1.y != p2.y)
    144.                 return p1.y < p2.y?1:0;
    145.         else
    146.                 return p1.x < p2.x?1:0;        
    147. }
    148. int cmp_yX (AcGePoint3d p1,AcGePoint3d p2)
    149. {
    150.         if (p1.y != p2.y)
    151.                 return p1.y < p2.y?1:0;
    152.         else
    153.                 return p1.x > p2.x?1:0;        
    154. }
    155. int cmp_yxz (AcGePoint3d p1,AcGePoint3d p2)
    156. {
    157.         if (p1.y != p2.y)
    158.                 return p1.y < p2.y?1:0;
    159.         else if (p1.x != p2.x)
    160.                 return p1.x < p2.x?1:0;
    161.         else
    162.                 return p1.z < p2.z?1:0;
    163. }
    164. int cmp_yXz (AcGePoint3d p1,AcGePoint3d p2)
    165. {
    166.         if (p1.y != p2.y)
    167.                 return p1.y < p2.y?1:0;
    168.         else if (p1.x != p2.x)
    169.                 return p1.x > p2.x?1:0;
    170.         else
    171.                 return p1.z < p2.z?1:0;
    172. }
    173. int cmp_yXZ (AcGePoint3d p1,AcGePoint3d p2)
    174. {
    175.         if (p1.y != p2.y)
    176.                 return p1.y < p2.y?1:0;
    177.         else if (p1.x != p2.x)
    178.                 return p1.x > p2.x?1:0;
    179.         else
    180.                 return p1.z > p2.z?1:0;
    181. }
    182. int cmp_yxZ (AcGePoint3d p1,AcGePoint3d p2)
    183. {
    184.         if (p1.y != p2.y)
    185.                 return p1.y < p2.y?1:0;
    186.         else if (p1.x != p2.x)
    187.                 return p1.x < p2.x?1:0;
    188.         else
    189.                 return p1.z > p2.z?1:0;
    190. }
    191. int cmp_YXZ (AcGePoint3d p1,AcGePoint3d p2)
    192. {
    193.         if (p1.y != p2.y)
    194.                 return p1.y > p2.y?1:0;
    195.         else if (p1.x != p2.x)
    196.                 return p1.x > p2.x?1:0;
    197.         else
    198.                 return p1.z > p2.z?1:0;
    199. }
    200. int cmp_YxZ (AcGePoint3d p1,AcGePoint3d p2)
    201. {
    202.         if (p1.y != p2.y)
    203.                 return p1.y > p2.y?1:0;
    204.         else if (p1.x != p2.x)
    205.                 return p1.x < p2.x?1:0;
    206.         else
    207.                 return p1.z > p2.z?1:0;
    208. }

    209. int cmp_Yxz (AcGePoint3d p1,AcGePoint3d p2)
    210. {
    211.         if (p1.y != p2.y)
    212.                 return p1.y > p2.y?1:0;
    213.         else if (p1.x != p2.x)
    214.                 return p1.x < p2.x?1:0;
    215.         else
    216.                 return p1.z < p2.z?1:0;
    217. }

    218. int cmp_YXz (AcGePoint3d p1,AcGePoint3d p2)
    219. {
    220.         if (p1.y != p2.y)
    221.                 return p1.y > p2.y?1:0;
    222.         else if (p1.x != p2.x)
    223.                 return p1.x > p2.x?1:0;
    224.         else
    225.                 return p1.z < p2.z?1:0;
    226. }


    227. int cmp_2d_xy (AcGePoint2d p1,AcGePoint2d p2)
    228. {
    229.         if (p1.x != p2.x)
    230.                 return p1.x < p2.x?1:0;
    231.         else
    232.                 return p1.y < p2.y?1:0;        
    233. }
    234. int cmp_2d_xY (AcGePoint2d p1,AcGePoint2d p2)
    235. {
    236.         if (p1.x != p2.x)
    237.                 return p1.x < p2.x?1:0;
    238.         else
    239.                 return p1.y > p2.y?1:0;        
    240. }
    241. int cmp_2d_Xy (AcGePoint2d p1,AcGePoint2d p2)
    242. {
    243.         if (p1.x != p2.x)
    244.                 return p1.x > p2.x?1:0;
    245.         else
    246.                 return p1.y < p2.y?1:0;        
    247. }
    248. int cmp_2d_XY (AcGePoint2d p1,AcGePoint2d p2)
    249. {
    250.         if (p1.x != p2.x)
    251.                 return p1.x > p2.x?1:0;
    252.         else
    253.                 return p1.y > p2.y?1:0;        
    254. }

    255. int cmp_2d_x (AcGePoint2d p1,AcGePoint2d p2)
    256. {
    257.         return p1.x < p2.x?1:0;
    258. }


    259. int cmp_2d_X(AcGePoint2d p1,AcGePoint2d p2)
    260. {
    261.         return p1.x > p2.x?1:0;
    262. }


    263. int cmp_2d_y (AcGePoint2d p1,AcGePoint2d p2)
    264. {
    265.         return p1.y < p2.y?1:0;
    266. }


    267. int cmp_2d_Y(AcGePoint2d p1,AcGePoint2d p2)
    268. {
    269.         return p1.y > p2.y?1:0;
    270. }

    271. int cmp_2d_YX(AcGePoint2d p1,AcGePoint2d p2)
    272. {
    273.         if (p1.y != p2.y)
    274.                 return p1.y > p2.y?1:0;
    275.         else
    276.                 return p1.x > p2.x?1:0;        
    277. }
    278. int cmp_2d_Yx (AcGePoint2d p1,AcGePoint2d p2)
    279. {
    280.         if (p1.y != p2.y)
    281.                 return p1.y > p2.y?1:0;
    282.         else
    283.                 return p1.x < p2.x?1:0;        
    284. }



    285. int cmp_2d_yx (AcGePoint2d p1,AcGePoint2d p2)
    286. {
    287.         if (p1.y != p2.y)
    288.                 return p1.y < p2.y?1:0;
    289.         else
    290.                 return p1.x < p2.x?1:0;        
    291. }



    292. int cmp_2d_yX (AcGePoint2d p1,AcGePoint2d p2)
    293. {
    294.         if (p1.y != p2.y)
    295.                 return p1.y < p2.y?1:0;
    296.         else
    297.                 return p1.x > p2.x?1:0;        
    298. }



    299. int cmp_up (int &a,int &b)
    300. {
    301.         return a<b?1:0;
    302. }



    303. int cmp_up (double &a,double &b)
    304. {
    305.         return a<b?1:0;
    306. }



    307. int cmp_down (const int &a, const int &b)
    308. {
    309.         return a>b?1:0;
    310. }


    311. int cmp_down (double &a,double &b)
    312. {
    313.         return a>b?1:0;
    314. }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-1-21 21:00 | 显示全部楼层
    1. 计算几何几何函数库
    2. -------------------------------------------------------------------------------------------------------------------------------
    3. 导引
    4. 1. 常量定义和包含文件
    5. 2. 基本数据结构
    6. 3. 精度控制
    7. ㈠ 点的基本运算
    8. 1. 平面上两点之间距离
    9. 2. 判断两点是否重合
    10. 3. 矢量叉乘
    11. 4. 矢量点乘
    12. 5. 判断点是否在线段上
    13. 6. 求一点饶某点旋转后的坐标
    14. 7. 求矢量夹角
    15. ㈡ 线段及直线的基本运算
    16. 1. 点与线段的关系
    17. 2. 求点到线段所在直线垂线的垂足
    18. 3. 点到线段的最近点
    19. 4. 点到线段所在直线的距离
    20. 5. 点到折线集的最近距离
    21. 6. 判断圆是否在多边形内
    22. 7. 求矢量夹角余弦
    23. 8. 求线段之间的夹角
    24. 9. 判断线段是否相交
    25. 10.判断线段是否相交但不交在端点处
    26. 11.求点关于某直线的对称点
    27. 12.判断两条直线是否相交及求直线交点
    28. 13.判断线段是否相交,如果相交返回交点
    29. ㈢ 多边形常用算法模块
    30. 1. 判断多边形是否简单多边形
    31. 2. 检查多边形顶点的凸凹性
    32. 3. 判断多边形是否凸多边形
    33. 4. 求多边形面积
    34. 5. 判断多边形顶点的排列方向
    35. 7. 射线法判断点是否在多边形内
    36. 8. 判断点是否在凸多边形内
    37. 9. 寻找点集的graham算法
    38. 10.寻找点集凸包的卷包裹法
    39. 11.凸包MelkMan算法的实现
    40. 12. 凸多边形的直径
    41. 13.求凸多边形的重心
    42. ===========================================================================
    43. 导引
    44. /* 需要包含的头文件 */
    45. #include <cmath >
    46. /* 常量定义 */
    47. const double INF = 1E200;
    48. const double EP = 1E-10;
    49. const int MAXV = 300;
    50. const double PI = 3.14159265;
    51. /* 基本几何结构 */
    52. struct POINT
    53. {
    54. double x;
    55. double y;
    56. POINT(double a=0, double b=0) { x=a; y=b;}
    57. };
    58. struct LINESEG
    59. {
    60. POINT s;
    61. POINT e;
    62. LINESEG(POINT a, POINT b) { s=a; e=b;}
    63. LINESEG() { }
    64. };
    65. // 直线的解析方程 a*x+b*y+c=0 为统一表示,约定 a>= 0
    66. struct LINE
    67. {
    68. double a;
    69. double b;
    70. double c;
    71. LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;}
    72. };
    73. //线段树
    74. struct LINETREE
    75. {
    76. }
    77. //浮点误差的处理
    78. int dblcmp(double d)
    79. {
    80. if(fabs(d)<EP)
    81. return 0 ;
    82. return (d>0) ?1 :-1 ;
    83. }
    84. <一>点的基本运算
    85. // 返回两点之间欧氏距离
    86. double dist(POINT p1,POINT p2)
    87. {
    88. return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );
    89. }
    90. // 判断两个点是否重合
    91. bool equal_point(POINT p1,POINT p2)
    92. {
    93. return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) );
    94. }
    95. /*(sp-op)*(ep-op)的叉积
    96. r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积
    97. r>0:sp在矢量op ep的顺时针方向;
    98. r=0:op sp ep三点共线;
    99. r<0: sp在矢量op ep的逆时针方向 */
    100. double multiply(POINT sp,POINT ep,POINT op)
    101. {
    102. return((sp.x-op.x)*(ep.y-op.y) - (ep.x-op.x)*(sp.y-op.y));
    103. }
    104. double amultiply(POINT sp,POINT ep,POINT op)
    105. {
    106. return fabs((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));
    107. }
    108. /*矢量(p1-op)和(p2-op)的点积
    109. r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积如果两个矢量都非零矢量
    110. r < 0: 两矢量夹角为锐角;
    111. r = 0:两矢量夹角为直角;
    112. r > 0: 两矢量夹角为钝角 */
    113. double dotmultiply(POINT p1,POINT p2,POINT p0)
    114. {
    115. return ((p1.x-p0.x)*(p2.x-p0.x) + (p1.y-p0.y)*(p2.y-p0.y));
    116. }
    117. /* 判断点p是否在线段l上
    118. 条件:(p在线段l所在的直线上)&& (点p在以线段l为对角线的矩形内) */
    119. bool online(LINESEG l,POINT p)
    120. {
    121. return ((multiply(l.e,p,l.s)==0)
    122. && ( ( (p.x-l.s.x) * (p.x-l.e.x) <=0 ) && ( (p.y-l.s.y)*(p.y-l.e.y) <=0 ) ) );
    123. }
    124. // 返回点p以点o为圆心逆时针旋转alpha(单位:弧度)后所在的位置
    125. POINT rotate(POINT o,double alpha,POINT p)
    126. {
    127. POINT tp;
    128. p.x -=o.x;
    129. p.y -=o.y;
    130. tp.x=p.x*cos(alpha) - p.y*sin(alpha)+o.x;
    131. tp.y=p.y*cos(alpha) + p.x*sin(alpha)+o.y;
    132. return tp;
    133. }
    134. /* 返回顶角在o点,起始边为os,终止边为oe的夹角(单位:弧度)
    135. 角度小于pi,返回正值
    136. 角度大于pi,返回负值
    137. 可以用于求线段之间的夹角 */
    138. double angle(POINT o,POINT s,POINT e)
    139. {
    140. double cosfi,fi,norm;
    141. double dsx = s.x - o.x;
    142. double dsy = s.y - o.y;
    143. double dex = e.x - o.x;
    144. double dey = e.y - o.y;
    145. cosfi=dsx*dex+dsy*dey;
    146. norm=(dsx*dsx+dey*dey)*(dex*dex+dey*dey);
    147. cosfi /= sqrt( norm );
    148. if (cosfi >= 1.0 ) return 0;
    149. if (cosfi <= -1.0 ) return -3.1415926;
    150. fi=acos(cosfi);
    151. if (dsx*dey-dsy*dex>0) return fi;// 说明矢量os 在矢量 oe的顺时针方向
    152. return -fi;
    153. }
    154. <二>线段及直线的基本运算
    155. /* 判断点C在线段AB所在的直线l上垂足P的与线段AB的关系
    156. 本函数是根据下面的公式写的,P是点C到线段AB所在直线的垂足
    157. AC dot AB
    158. r = ----------------------
    159. ||AB||^2
    160. (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    161. = ----------------------------------------------------
    162. L^2
    163. r has the following meaning:
    164. r=0 P = A
    165. r=1 P = B
    166. r<0 P is on the backward extension of AB
    167. r>1 P is on the forward extension of AB
    168. 0<r<1 P is interior to AB
    169. */
    170. double relation(POINT c,LINESEG l)
    171. {
    172. LINESEG tl;
    173. tl.s=l.s;
    174. tl.e=c;
    175. return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e));
    176. }
    177. // 求点C到线段AB所在直线的垂足 P
    178. POINT perpendicular(POINT p,LINESEG l)
    179. {
    180. double r=relation(p,l);
    181. POINT tp;
    182. tp.x=l.s.x+r*(l.e.x-l.s.x);
    183. tp.y=l.s.y+r*(l.e.y-l.s.y);
    184. return tp;
    185. }
    186. /* 求点p到线段l的最短距离
    187. 返回线段上距该点最近的点np 注意:np是线段l上到点p最近的点,不一定是垂足 */
    188. double ptolinesegdist(POINT p,LINESEG l,POINT &np)
    189. {
    190. double r=relation(p,l);
    191. if(r<0)
    192. {
    193. np=l.s;
    194. return dist(p,l.s);
    195. }
    196. if(r>1)
    197. {
    198. np=l.e;
    199. return dist(p,l.e);
    200. }
    201. np=perpendicular(p,l);
    202. return dist(p,np);
    203. }
    204. // 求点p到线段l所在直线的距离
    205. //请注意本函数与上个函数的区别
    206. double ptoldist(POINT p,LINESEG l)
    207. {
    208. return abs(multiply(p,l.e,l.s))/dist(l.s,l.e);
    209. }
    210. /* 计算点到折线集的最近距离,并返回最近点.
    211. 注意:调用的是ptolineseg()函数 */
    212. double ptopointset(int vcount, POINT pointset[], POINT p, POINT &q)
    213. {
    214. int i;
    215. double cd=double(INF),td;
    216. LINESEG l;
    217. POINT tq,cq;
    218. for(i=0;i<vcount-1;i++)
    219. {
    220. l.s=pointset[i];
    221. l.e=pointset[i+1];
    222. td=ptolinesegdist(p,l,tq);
    223. if(td<cd)
    224. {
    225. cd=td;
    226. cq=tq;
    227. }
    228. }
    229. q=cq;
    230. return cd;
    231. }
    232. /* 判断圆是否在多边形内*/
    233. bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[])
    234. {
    235. POINT q;
    236. double d;
    237. q.x=0;
    238. q.y=0;
    239. d=ptopointset(vcount,polygon,center,q);
    240. if(d<radius||fabs(d-radius)<EP) return true;
    241. else return false;
    242. }
    243. /* 返回两个矢量l1和l2的夹角的余弦 (-1 ~ 1)
    244. 注意:如果想从余弦求夹角的话,注意反余弦函数的值域是从 0到pi */
    245. double cosine(LINESEG l1,LINESEG l2)
    246. {
    247. return(((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x)+(l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) );
    248. }
    249. // 返回线段l1与l2之间的夹角
    250. //单位:弧度 范围(-pi,pi)
    251. double lsangle(LINESEG l1,LINESEG l2)
    252. {
    253. POINT o,s,e;
    254. o.x=o.y=0;
    255. s.x=l1.e.x-l1.s.x;
    256. s.y=l1.e.y-l1.s.y;
    257. e.x=l2.e.x-l2.s.x;
    258. e.y=l2.e.y-l2.s.y;
    259. return angle(o,s,e);
    260. }
    261. //判断线段u和v相交(包括相交在端点处)
    262. bool intersect(LINESEG u,LINESEG v)
    263. {
    264. return ( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&& //排斥实验
    265. (max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&&
    266. (max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&&
    267. (max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&&
    268. (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&& //跨立实验
    269. (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0));
    270. }
    271. // 判断线段u和v相交(不包括双方的端点)
    272. bool intersect_A(LINESEG u,LINESEG v)
    273. {
    274. return ((intersect(u,v)) &&
    275. (!online(u,v.s)) &&
    276. (!online(u,v.e)) &&
    277. (!online(v,u.e)) &&
    278. (!online(v,u.s)));
    279. }
    280. // 判断线段v所在直线与线段u相交
    281. 方法:判断线段u是否跨立线段v
    282. bool intersect_l(LINESEG u,LINESEG v)
    283. {
    284. return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0;
    285. }
    286. // 根据已知两点坐标,求过这两点的直线解析方程: a*x+b*y+c = 0 (a >= 0)
    287. LINE makeline(POINT p1,POINT p2)
    288. {
    289. LINE tl;
    290. int sign = 1;
    291. tl.a=p2.y-p1.y;
    292. if(tl.a<0)
    293. {
    294. sign = -1;
    295. tl.a=sign*tl.a;
    296. }
    297. tl.b=sign*(p1.x-p2.x);
    298. tl.c=sign*(p1.y*p2.x-p1.x*p2.y);
    299. return tl;
    300. }
    301. // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200
    302. double slope(LINE l)
    303. {
    304. if(abs(l.a) < 1e-20)return 0;
    305. if(abs(l.b) < 1e-20)return INF;
    306. return -(l.a/l.b);
    307. }
    308. // 返回直线的倾斜角alpha ( 0 - pi)
    309. // 注意:atan()返回的是 -PI/2 ~ PI/2
    310. double alpha(LINE l)
    311. {
    312. if(abs(l.a)< EP)return 0;
    313. if(abs(l.b)< EP)return PI/2;
    314. double k=slope(l);
    315. if(k>0)
    316. return atan(k);
    317. else
    318. return PI+atan(k);
    319. }
    320. // 求点p关于直线l的对称点
    321. POINT symmetry(LINE l,POINT p)
    322. {
    323. POINT tp;
    324. tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
    325. tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
    326. return tp;
    327. }
    328. // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交,返回true,且返回交点p
    329. bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1,L2
    330. {
    331. double d=l1.a*l2.b-l2.a*l1.b;
    332. if(abs(d)<EP) // 不相交
    333. return false;
    334. p.x = (l2.c*l1.b-l1.c*l2.b)/d;
    335. p.y = (l2.a*l1.c-l1.a*l2.c)/d;
    336. return true;
    337. }
    338. // 如果线段l1和l2相交,返回true且交点由(inter)返回,否则返回false
    339. bool intersection(LINESEG l1,LINESEG l2,POINT &inter)
    340. {
    341. LINE ll1,ll2;
    342. ll1=makeline(l1.s,l1.e);
    343. ll2=makeline(l2.s,l2.e);
    344. if(lineintersect(ll1,ll2,inter)) return online(l1,inter);
    345. else return false;
    346. }
    347. <三> 多边形常用算法模块
    348. 如果无特别说明,输入多边形顶点要求按逆时针排列
    349. // 返回多边形面积(signed);
    350. // 输入顶点按逆时针排列时,返回正值;否则返回负值
    351. double area_of_polygon(int vcount,POINT polygon[])
    352. {
    353. int i;
    354. double s;
    355. if (vcount<3)
    356. return 0;
    357. s=polygon[0].y*(polygon[vcount-1].x-polygon[1].x);
    358. for (i=1;i<vcount;i++)
    359. s+=polygon[i].y*(polygon[(i-1)].x-polygon[(i+1)%vcount].x);
    360. return s/2;
    361. }
    362. // 判断顶点是否按逆时针排列
    363. // 如果输入顶点按逆时针排列,返回true
    364. bool isconterclock(int vcount,POINT polygon[])
    365. {
    366. return area_of_polygon(vcount,polygon)>0;
    367. }
    368. /*射线法判断点q与多边形polygon的位置关系
    369. 要求polygon为简单多边形,顶点时针排列
    370. 如果点在多边形内: 返回0
    371. 如果点在多边形边上:返回1
    372. 如果点在多边形外: 返回2 */
    373. int insidepolygon(POINT q)
    374. {
    375. int c=0,i,n;
    376. LINESEG l1,l2;
    377. l1.s=q; l1.e=q;l1.e.x=double(INF);
    378. n=vcount;
    379. for (i=0;i<vcount;i++)
    380. {
    381. l2.s=Polygon[i];
    382. l2.e=Polygon[(i+1)%vcount];
    383. double ee= Polygon[(i+2)%vcount].x;
    384. double ss= Polygon[(i+3)%vcount].y;
    385. if(online(l2,q))
    386. return 1;
    387. if(intersect_A(l1,l2))
    388. c++; // 相交且不在端点
    389. if(online(l1,l2.e)&& !online(l1,l2.s) && l2.e.y>l2.e.y)
    390. c++;//l2的一个端点在l1上且该端点是两端点中纵坐标较大的那个
    391. if(!online(l1,l2.e)&& online(l1,l2.s) && l2.e.y<l2.e.y)
    392. c++;//忽略平行边
    393. }
    394. if(c%2 == 1)
    395. return 0;
    396. else
    397. return 2;
    398. }
    399. //判断点q在凸多边形polygon内
    400. // 点q是凸多边形polygon内[包括边上]时,返回true
    401. // 注意:多边形polygon一定要是凸多边形
    402. bool InsideConvexPolygon(int vcount,POINT polygon[],POINT q)
    403. {
    404. POINT p;
    405. LINESEG l;
    406. int i;
    407. p.x=0; p.y=0;
    408. for(i=0;i<vcount;i++) // 寻找一个肯定在多边形polygon内的点p:多边形顶点平均值
    409. {
    410. p.x+=polygon[i].x;
    411. p.y+=polygon[i].y;
    412. }
    413. p.x /= vcount;
    414. p.y /= vcount;
    415. for(i=0;i<vcount;i++)
    416. {
    417. l.s=polygon[i];
    418. l.e=polygon[(i+1)%vcount];
    419. if(multiply(p,l.e,l.s)*multiply(q,l.e,l.s)<0)
    420. /* 点p和点q在边l的两侧,说明点q肯定在多边形外 */
    421. return false;
    422. }
    423. return true;
    424. }
    425. /*寻找凸包的graham 扫描法
    426. PointSet为输入的点集;
    427. ch为输出的凸包上的点集,按照逆时针方向排列;
    428. n为PointSet中的点的数目
    429. len为输出的凸包上的点的个数 */
    430. void Graham_scan(POINT PointSet[],POINT ch[],int n,int &len)
    431. {
    432. int i,j,k=0,top=2;
    433. POINT tmp;
    434. // 选取PointSet中y坐标最小的点PointSet[k],如果这样的点有多个,则取最左边的一个
    435. for(i=1;i<n;i++)
    436. if ( PointSet[i].y<PointSet[k].y || (PointSet[i].y==PointSet[k].y)
    437. && (PointSet[i].x<PointSet[k].x) )
    438. k=i;
    439. tmp=PointSet[0];
    440. PointSet[0]=PointSet[k];
    441. PointSet[k]=tmp; // 现在PointSet中y坐标最小的点在PointSet[0]
    442. for (i=1;i<n-1;i++) /* 对顶点按照相对PointSet[0]的极角从小到大进行排序,极角相同
    443. 的按照距离PointSet[0]从近到远进行排序 */
    444. {
    445. k=i;
    446. for (j=i+1;j<n;j++)
    447. if ( multiply(PointSet[j],PointSet[k],PointSet[0])>0 || // 极角更小
    448. (multiply(PointSet[j],PointSet[k],PointSet[0])==0) && /*极角相等,距离更短 */ dist(PointSet[0],PointSet[j])<dist(PointSet[0],PointSet[k]) )
    449. k=j;
    450. tmp=PointSet[i];
    451. PointSet[i]=PointSet[k];
    452. PointSet[k]=tmp;
    453. }
    454. ch[0]=PointSet[0];
    455. ch[1]=PointSet[1];
    456. ch[2]=PointSet[2];
    457. for (i=3;i<n;i++)
    458. {
    459. while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;
    460. ch[++top]=PointSet[i];
    461. }
    462. len=top+1;
    463. }
    464. // 卷包裹法求点集凸壳,参数说明同graham算法
    465. void ConvexClosure(POINT PointSet[],POINT ch[],int n,int &len)
    466. {
    467. int top=0,i,index,first;
    468. double curmax,curcos,curdis;
    469. POINT tmp;
    470. LINESEG l1,l2;
    471. bool use[MAXV];
    472. tmp=PointSet[0];
    473. index=0;
    474. // 选取y最小点,如果多于一个,则选取最左点
    475. for(i=1;i<n;i++)
    476. {
    477. if(PointSet[i].y<tmp.y||PointSet[i].y == tmp.y&&PointSet[i].x<tmp.x)
    478. {
    479. index=i;
    480. }
    481. use[i]=false;
    482. }
    483. tmp=PointSet[index];
    484. first=index;
    485. use[index]=true;
    486. index=-1;
    487. ch[top++]=tmp;
    488. tmp.x-=100;
    489. l1.s=tmp;
    490. l1.e=ch[0];
    491. l2.s=ch[0];
    492. while(index!=first)
    493. {
    494. curmax=-100;
    495. curdis=0;
    496. // 选取与最后一条确定边夹角最小的点,即余弦值最大者
    497. for(i=0;i<n;i++)
    498. {
    499. if(use[i])continue;
    500. l2.e=PointSet[i];
    501. curcos=cosine(l1,l2); // 根据cos值求夹角余弦,范围在 (-1 -- 1 )
    502. if(curcos>curmax || fabs(curcos-curmax)<1e-6 && dist(l2.s,l2.e)>curdis)
    503. {
    504. curmax=curcos;
    505. index=i;
    506. curdis=dist(l2.s,l2.e);
    507. }
    508. }
    509. use[first]=false; //清空第first个顶点标志,使最后能形成封闭的hull
    510. use[index]=true;
    511. ch[top++]=PointSet[index];
    512. l1.s=ch[top-2];
    513. l1.e=ch[top-1];
    514. l2.s=ch[top-1];
    515. }
    516. len=top-1;
    517. }
    518. // 求凸多边形的重心,要求输入多边形按逆时针排序
    519. POINT gravitycenter(int vcount,POINT polygon[])
    520. {
    521. POINT tp;
    522. double x,y,s,x0,y0,cs,k;
    523. x=0;y=0;s=0;
    524. for(int i=1;i<vcount-1;i++)
    525. {
    526. x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3;
    527. y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求当前三角形的重心
    528. cs=multiply(polygon[i],polygon[i+1],polygon[0])/2;
    529. //三角形面积可以直接利用该公式求解
    530. if(abs(s)<1e-20)
    531. {
    532. x=x0;y=y0;s+=cs;continue;
    533. }
    534. k=cs/s; //求面积比例
    535. x=(x+k*x0)/(1+k);
    536. y=(y+k*y0)/(1+k);
    537. s += cs;
    538. }
    539. tp.x=x;
    540. tp.y=y;
    541. return tp;
    542. }
    543. /*所谓凸多边形的直径,即凸多边形任两个顶点的最大距离。下面的算法
    544. 仅耗时O(n),是一个优秀的算法。 输入必须是一个凸多边形,且顶点
    545. 必须按顺序(顺时针、逆时针均可)依次输入。若输入不是凸多边形
    546. 而是一般点集,则要先求其凸包。 就是先求出所有跖对,然后求出每
    547. 个跖对的距离,取最大者。点数要多于5个*/
    548. void Diameter(POINT ch[],int n,double &dia)
    549. {
    550. int znum=0,i,j,k=1;
    551. int zd[MAXV][2];
    552. double tmp;
    553. while(amultiply(ch[0],ch[k+1],ch[n-1]) > amultiply(ch[0],ch[k],ch[n-1])-EP)
    554. k++;
    555. i=0;
    556. j=k;
    557. while(i<=k && j<n)
    558. {
    559. zd[znum][0]=i;
    560. zd[znum++][1]=j;
    561. while(amultiply(ch[i+1],ch[j+1],ch[i]) > amultiply(ch[i+1],ch[j],ch[i]) – EP
    562. && j< n-1)
    563. {
    564. zd[znum][0]=i;
    565. zd[znum++][1]=j;
    566. j++;
    567. }
    568. i++;
    569. }
    570. dia=-1.0;
    571. for(i=0;i<znum;i++)
    572. {
    573. printf("%d %d/n",zd[i][0],zd[i][1]);
    574. tmp=dist(ch[zd[i][0]],ch[zd[i][1]]);
    575. if(dia<tmp)
    576. dia=tmp;
    577. }
    578. }
    复制代码
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    关闭

    推荐膜材品牌上一条 /6 下一条

    进口膜材 国产膜材 pvdf膜材ptfe膜材ETFE膜材
    最好的膜结构公司 一级膜结构资质 膜结构一级资质
    膜结构设计-膜结构十大品牌-etfe设计-充气膜结构
    诺科膜结构
    遨都膜结构设计
    中国膜结构网
    中国空间膜结构

    QQ|申请友链|手机版|中国膜结构论坛