兩條相交的線段AI,BI，交點為I
圓O半徑為R，與兩線相切，切點為M,N，M,N分別在AI,BI上。
點C是M,N兩點之間的劣弧的中點
現在已知: A(x0,y0)，B(x1,y1)，C(x2,y2)三點的坐標，求圓心O 的坐標xy
目前我自己寫了一個演算法，是按微積分的方式求得的解，感覺應該是有決議解的，請CSDN的大神們幫忙看看
uj5u.com熱心網友回復：
我的微積分演算法代碼，只作測驗使用，不完善之處，敬請指正。
// pppl.cpp : 定義控制臺應用程式的入口點。
//
/*********************************************************************************
* Copyright(C),;
* FileName: pppl.cpp
* Author: lizhq;
* Date: 2019-1-18
* Description: 已知相交線段的兩個發散側端點及交點側倒角圓弧的中點及圓弧半徑，
求該曲線(兩條線段+圓弧)，
轉化為求圓心坐標，圓心坐標求出后，求曲線就很簡單了。
*********************************************************************************/
#include "math.h"
#include <vector>
#define min(a,b) (((a) < (b)) ? (a) : (b))
#define max(a,b) (((a) > (b)) ? (a) : (b))
#define eps (1e-12)
#define PI (3.1415926535897932384626)
#define R2D(rad) (rad / PI * 180.0)
#define D2R(ang) (ang / 180.0 * PI)
class CLPoit
{
public:
double x;
double y;
CLPoit()
{
y = 0.;
x = 0.;
}
CLPoit(double xx, double yy)
{
y = yy;
x = xx;
}
};
class CLArc
{
public:
double x;
double y;
double r;
double staRadian;
double sweepRadian;
CLArc()
{
y = 0.;
x = 0.;
r = 0.;
staRadian = 0.;
sweepRadian = 0.;
}
CLArc(const CLPoit& c, double rr, double sta, double swp)
{
x = c.x;
y = c.y;
r = rr;
staRadian = sta;
sweepRadian = swp;
}
CLArc(double xx, double yy, double rr, double sta, double swp)
{
x = xx;
y = yy;
r = rr;
staRadian = sta;
sweepRadian = swp;
}
void set(const CLPoit& c, double rr, double sta, double swp)
{
x = c.x;
y = c.y;
r = rr;
staRadian = sta;
sweepRadian = swp;
}
void set(double xx, double yy, double rr, double sta, double swp)
{
x = xx;
y = yy;
r = rr;
staRadian = sta;
sweepRadian = swp;
}
void getLinePts(std::vector<CLPoit>& pts, int count)
{
pts.clear();
if (count < 2)
return;
CLPoit pt;
double argRad = sweepRadian / (count - 1);
for (int i = 0; i < count; i++)
{
pt.x = r*cos((staRadian + argRad*double(i))) + x;
pt.y = r*sin((staRadian + argRad*double(i))) + y;
pts.push_back(pt);
}
}
};
double dis(const CLPoit& a, const CLPoit& b)
{
return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
bool getQieDian(CLPoit& tangPt, const CLPoit& p, const CLPoit& o, double r, bool nsz)
{
double m = p.x, n = p.y;
double a = o.x, b = o.y;
// 點到圓心距離的平方
double d2 = (m - a) * (m - a) + (n - b) * (n - b);
// 點到圓心距離
double d = sqrt(d2);
if (d < r - eps)
return false;
if (d < r + eps)
{
tangPt = p;
return true;
}
// 半徑的平方
double r2 = r * r;
// 點到切點距離
double l = sqrt(d2 - r2);
// 點->圓心的單位向量
double x0 = (a - m) / d;
double y0 = (b - n) / d;
// 計算切線與點心連線的夾角
double f = asin(r / d);
// 向正反兩個方向旋轉單位向量
if (!nsz) f = -f;
double x1 = x0 * cos(f) - y0 * sin(f);
double y1 = x0 * sin(f) + y0 * cos(f);
// 得到新座標
tangPt.x = x1 * l + m;
tangPt.y = y1 * l + n;
return true;
}
double getRadian(const CLPoit& p1, const CLPoit& p2)
{
double radian_temp;
double xx, yy;
xx = p2.x - p1.x;
yy = p2.y - p1.y;
if (xx == 0.0)
radian_temp = PI / 2.0;
else
radian_temp = atan(fabs(yy / xx));
if ((xx < 0.0) && (yy >= 0.0))
radian_temp = PI - radian_temp;
else if ((xx < 0.0) && (yy < 0.0))
radian_temp = PI + radian_temp;
else if ((xx >= 0.0) && (yy < 0.0))
radian_temp = PI * 2.0 - radian_temp;
return (radian_temp);
}
int getMinValIdx(double& minVal, const std::vector<CLPoit>& pts, const CLPoit& ptL, const CLPoit& ptM, const CLPoit& ptR, double r, bool nszxz)
{
int idx = -1;
minVal = DBL_MAX;
for (size_t i = 0; i < pts.size(); i++)
{
auto& ori = pts[i];
CLPoit qd1, qd2;
if (!getQieDian(qd1, ptL, ori, r, nszxz) || !getQieDian(qd2, ptR, ori, r, !nszxz))
continue;
double l1 = dis(qd1, ptM);
double l2 = dis(qd2, ptM);
double dLen = fabs(l1 - l2);
if (dLen < minVal)
{
minVal = dLen;
idx = i;
}
}
return idx;
}
void getAllLen(CLPoit& ori, const CLPoit& ptL, const CLPoit& ptM, const CLPoit& ptR, double r)
{
std::vector<double> lens;
double staRadian = getRadian(ptM, ptL);
double swRadian = getRadian(ptM, ptR) - staRadian;
if (swRadian > PI) swRadian -= 2 * PI;
else if (swRadian < -PI) swRadian += 2 * PI;
if (fabs(swRadian) < eps || fabs(fabs(swRadian) - PI) < eps)
{
//三點共線，無需找圓心
ori = ptM;
return;
}
bool nszxz = (swRadian > 0);
CLArc arc(ptM, r, staRadian, swRadian);
double deg = fabs(R2D(swRadian));
int count = (deg > 10.0) ? int(deg + 0.5) : 10;
std::vector<CLPoit> pts;
arc.getLinePts(pts, count);
for (int i = 0; i < count; i++)
{
ori = pts[i];
double dLen = -1.0;
CLPoit qd1, qd2;
if( getQieDian(qd1, ptL, ori, r, nszxz) && getQieDian(qd2, ptR, ori, r, !nszxz))
{
double l1 = dis(qd1, ptM);
double l2 = dis(qd2, ptM);
dLen = fabs(l1 - l2);
}
lens.push_back(dLen);
}
}
// 使用如下測驗資料測驗時出現死回圈
// double r = 8000;
// CLPoit p1(-48248.5749, -14687.7393);
// CLPoit p2(-77030.8830, -6882.3034);
// CLPoit p3(-78825.9401, 6136.8982);
// 經分析是swRadian值過小時（e-15數量級），迭代程序中其值不再減小，造成死回圈
// 推斷可能與編譯環境及作業系統位數相關，在不清楚具體原因的情況下，為避免死回圈
// 應當增加判斷，若swRadian的絕對值不再繼續減小時，應當跳出while回圈
int getOriPoint(CLPoit& ori, const CLPoit& ptL, const CLPoit& ptM, const CLPoit& ptR, double r)
{
//getAllLen(ori, ptL, ptM, ptR, r);
ori = ptM;
double staRadian = getRadian(ptM, ptL);
double swRadian = getRadian(ptM, ptR) - staRadian;
if (swRadian > PI) swRadian -= 2 * PI;
else if (swRadian < -PI) swRadian += 2 * PI;
if (fabs(swRadian) < eps || fabs(fabs(swRadian) - PI) < eps)
return 1; //三點共線，無需找圓心
// 點ptL到最終圓弧的切線的旋轉方向
bool nszxz = (swRadian > 0);
//將圓弧劃分成最少10等份（11個點）;
std::vector<CLPoit> pts;
CLArc arc(ptM, r, staRadian, swRadian);
double deg = fabs(R2D(swRadian));
int count = (deg > 10.0) ? int(deg + 0.5) + 1 : 11;
arc.getLinePts(pts, count);
if (pts.empty())
return 3;
double lastSWRadian = swRadian;
double minVal = 0.0;
int idx, perIdx, nexIdx;
while (dis(pts.front(), pts.back()) > eps)
{
// 獲取最小值及對應點的索引值;
idx = getMinValIdx(minVal, pts, ptL, ptM, ptR, r, nszxz);
if (idx < 0)
return 2;
ori = pts[idx];
if (minVal < eps)
{
break;
}
//取最小值兩側的點，將圓弧劃分成4等份(5個點)進行迭代計算
perIdx = max(idx - 1, 0);
nexIdx = min(idx + 1, pts.size() - 1);
staRadian = getRadian(ptM, pts[perIdx]);
swRadian = getRadian(ptM, pts[nexIdx]) - staRadian;
if (swRadian > PI) swRadian -= 2 * PI;
else if (swRadian < -PI) swRadian += 2 * PI;
// 精度問題造成回圈時，跳出回圈(此時不判斷精度，直接判==);
if (lastSWRadian == swRadian)
break;
lastSWRadian = swRadian;
arc.set(ptM, r, staRadian, swRadian);
arc.getLinePts(pts, 5);
}
//校驗;
CLPoit qd1, qd2;
if (!getQieDian(qd1, ptL, ori, r, nszxz) || !getQieDian(qd2, ptR, ori, r, !nszxz))
return 2;
minVal = fabs(dis(qd1, ptM) - dis(qd2, ptM));
if (minVal < 1e-6)//校驗精度值可降低要求;
{
return 0;
}
return 2;
}
void Test()
{
double r = 8000;
CLPoit p1(-48248.5749, -14687.7393);
CLPoit p2(-77030.8830, -6882.3034);
CLPoit p3(-78825.9401, 6136.8982);
CLPoit ori;
int retVal = getOriPoint(ori, p1, p2, p3, r);
}
int main()
{
Test();
return 0;
}
uj5u.com熱心網友回復：
你打A點 B點什么意思uj5u.com熱心網友回復：
這個先設點建立方程組，表示出圓心坐標，得到橫縱坐標的運算式，帶入上述點的坐標值，就出來了uj5u.com熱心網友回復：
你是想提問還是干神魔，代碼怎么會那摩多uj5u.com熱心網友回復：
缺條件應該無解才對呀親。過O點隨便做一個圓，只要確保A, B點在圓外，則經過A, B點做圓的切線一定滿足題目的要求，可是根據作圖的方法，真正在畫圓的時候只用到了C的坐標，圓心O的位置完全可以隨便定。是我漏掉了什么嗎？
uj5u.com熱心網友回復：
半徑是定的，而且C是兩個相切點對應的弧線的中點，所以O不是可以隨便定的
uj5u.com熱心網友回復：
O的x，y作為未知數：用x，y以及A，B的坐標得到切線AM，BN的方程（幾何知識忘光了，應該可以得到吧？）
得到OC的方程
通過這三條直線可以兩兩算交點，有3種方程計算點I的運算式，這三種運算式肯定得相等，這樣得到2個方程式，而未知數只有2個，因該可以求解x和y。
uj5u.com熱心網友回復：
按說應該是可以求出決議解的，我也是跟你一樣，幾何知識忘光了，剛開始的時候就跟你想的一樣，算了半天，帶來帶去搞暈了，所以才發出來。
我就是lz，我lz這個號因發了一次主機屋的廣告，被禁言了，一直不得解，坑爹的主機屋，大家引以為戒吧。
uj5u.com熱心網友回復：
6樓正解，我的失誤，忘記說半徑是已知的定值。因C是弧MN 的中點，所以圓心應該只有一個值
uj5u.com熱心網友回復：
代碼是用微積分的方法求得的解，發出來希望得到決議解；
uj5u.com熱心網友回復：
今天有點時間大概算了下，如果樓主可以就自行化解，老夫過幾日再看看
uj5u.com熱心網友回復：
https://blog.csdn.net/bincam/article/details/109436568請各位大哥指點，多謝！
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/198806.html
標籤：數據結構與算法
上一篇：求助cmake編譯
下一篇：逆向之制作掃雷外掛——003