求过圆心直线与圆的两个交点

来源:互联网 发布:温职院网络课程 编辑:程序博客网 时间:2024/04/18 21:47

主要是注意所使用的数据类型。

之前用的是float,出现了一些意外,而且花费了我不少时间来反复验证、推导,

做了很多的无用功,而且,反复推导得出来的计算步骤并没有什么不牢靠的地方。

然后计算得到的结果却是让人如此之不省心,梗的我闷得慌。

今天上午发来了一贴,多位朋友各抒己见,

总算是让我发现了一些不足的地方,首当其冲的是一个变量弄错了,

导致大批的计算失准。

后来修正了这个bug以后,还是会出现计算不精确的地方。

再后来便将涉及的所有成员变量由float 纠正为 double 类型,

计算精度果然得到了提高,失准的地方再次被干掉。

这次给自己的教训就是:

涉及到精度比较高的数值运算的时候,还是得统统用 double。

之前还以为 float 已经比较不错,能够满足基本的需求了,

经过这次我总算是懂了,double的存在离我并不遥远。

这个问题堵了我比较久了,大概也有快10个月了,当时没解决就规避之没去用了,

今天能够解决这个遗留已久的问题,真是让人心情愉快!

下面贴出 Objective-C 和 Java 的相关代码:

Objective-C 部分(核心代码摘录)

[cpp] view plaincopy
  1. /** 已知两点,求过该两点的直线表达式~ */  
  2. - (BYLine) getLine:(b2Vec2)p1 anotherPoint:(b2Vec2)p2 {  
  3.     BYLine line;  
  4.     if((p1.x - p2.x) != 0) {  
  5.         line.kExists = true;  
  6.         line.k = (p1.y - p2.y) / (p1.x - p2.x);  
  7.         line.b = p1.y - line.k * p1.x;  
  8.     } else {  
  9.         line.kExists = false;  
  10.         line.extraX = p1.x;  
  11.     }  
  12.     return line;  
  13. }  
  14. /** 已知一点和直线斜率,求该直线的表达式~ */  
  15. - (BYLine) getLine:(b2Vec2)point kParam:(double)kParam {  
  16.     BYLine line;  
  17.     line.kExists = true;  
  18.     line.k = kParam;  
  19.     line.b = point.y - kParam * point.x;  
  20.     return line;  
  21. }  
  22. - (double) getDistanceBetween2Points:(b2Vec2)p0 anotherPoint:(b2Vec2)p1 {  
  23.     return sqrt(pow(p0.y - p1.y, 2) + pow(p0.x - p1.x, 2));  
  24. }  
  25. /** 获取一条直线上距离某点一定距离的两个点~ */  
  26. - (b2Vec2*) get2Points:(BYLine)ln p:(b2Vec2)point pw:(double)pathWidth {  
  27.     b2Vec2* target = new b2Vec2[2];  
  28.     double circleRadius = pathWidth / 2;  
  29.       
  30.     if(ln.k != 0) {  
  31.         // 斜率存在且不为 0~  
  32.         double kOfNewLine = -1 / ln.k;  
  33.         BYLine newLine = [self getLine:point kParam:kOfNewLine];  
  34.           
  35.         // 经过数学运算,得出二元一次方程组的表达式  
  36.         double A = pow(newLine.k, 2) + 1;  
  37.         double B = 2 * (newLine.k * newLine.b - newLine.k * point.y - point.x);  
  38.         double C = pow(point.x, 2) + pow((newLine.b - point.y), 2) - pow(circleRadius, 2);  
  39.         double delta = pow(B, 2) - 4 * A * C;  
  40.           
  41.         if(delta < 0) {    // 经实践检验有一定几率走入该分支,必须做特殊化处理~  
  42.             NSLog(@"竟然会无解,他妈的怎么回事儿啊!");  
  43.             target[0] = b2Vec2(point.x, point.y - circleRadius);  
  44.             target[1] = b2Vec2(point.x, point.y + circleRadius);  
  45.         } else {  
  46.             double x1 = (-B + sqrt(delta)) / (2 * A);  
  47.             double y1 = newLine.k * x1 + newLine.b;  
  48.             target[0] = b2Vec2(x1, y1);  
  49.               
  50.             double x2 = (-B - sqrt(delta)) / (2 * A);  
  51.             double y2 = newLine.k * x2 + newLine.b;  
  52.             target[1] = b2Vec2(x2, y2);  
  53.         }  
  54.     } else {  
  55.         // 斜率存在且为 0~  
  56.         target[0] = b2Vec2(point.x, point.y - circleRadius);  
  57.         target[1] = b2Vec2(point.x, point.y + circleRadius);  
  58.     }  
  59.     NSLog(@"离中心点的距离为:%f", [self getDistanceBetween2Points:target[0] anotherPoint:point]);  
  60.     return target;  
  61. }  
  62. // 绘制触摸点到移动点的轨迹,1个像素~  
  63. - (void) drawTouchPath {  
  64.     if(_mouseDown) {  
  65.         // 已知(2等分,用分数表示~)  
  66.         b2Vec2 pStart = _touchSegment.p1;  
  67.         b2Vec2 pEnd = _touchSegment.p2;  
  68.           
  69.         // 推出  
  70.         b2Vec2 pMiddle = b2Vec2((pStart.x + pEnd.x) / 2, (pStart.y + pEnd.y) / 2);  
  71.         float pathLength = [self getDistanceBetween2Points:pStart anotherPoint:pEnd];  
  72.           
  73.         // 设置触摸轨迹的宽度~  
  74.         float pathWidth = pathLength / 3.0f;  
  75.         if(pathWidth > TOUCH_PATH_MAX_WIDTH) {  
  76.             pathWidth = TOUCH_PATH_MAX_WIDTH;  
  77.         }  
  78.           
  79.         b2Vec2* result;  
  80.         BYLine expFunc = [self getLine:pStart anotherPoint:pEnd];  
  81.         if(expFunc.kExists) {   // 斜率存在~  
  82.             result = [self get2Points:expFunc p:pMiddle pw:pathWidth];  
  83.         } else {                // 斜率不存在~  
  84.             result = new b2Vec2[2];  
  85.             result[0] = b2Vec2(pMiddle.x - pathWidth / 2, pMiddle.y);  
  86.             result[1] = b2Vec2(pMiddle.x + pathWidth / 2, pMiddle.y);  
  87.         }  
  88.           
  89.         b2Vec2 finalResult[5];  
  90.         finalResult[0] = pStart;  
  91.         finalResult[1] = result[0];  
  92.         finalResult[2] = pEnd;  
  93.         finalResult[3] = result[1];  
  94.         finalResult[4] = pStart;  
  95.   
  96.         // 绘制白色内容物~  
  97.         glColor4f(1.0f, 1.0f, 1.0f, 1.0f);  
  98.         glVertexPointer(2, GL_FLOAT, 0, finalResult);  
  99.         glDrawArrays(GL_TRIANGLE_STRIP, 0, 5);  
  100.     }  
  101. }  

Java 部分(部件齐全,能直接拿来跑的)

[java] view plaincopy
  1. package org.bruce.vertices.controller.geometry;  
  2.   
  3. /** 
  4.  * @author BruceYang 
  5.  * 对点的抽象~ 
  6.  */  
  7. public class CGPoint {  
  8.     public double x;  
  9.     public double y;  
  10.       
  11.     public CGPoint() {  
  12.           
  13.     }  
  14.     public CGPoint(double x, double y) {  
  15.         this.x = x;  
  16.         this.y = y;  
  17.     }  
  18.       
  19.     @Override  
  20.     public String toString() {  
  21.         return "x=" + this.x + ", y=" + this.y;  
  22.     }  
  23. }  
*****************************************************

[java] view plaincopy
  1. package org.bruce.vertices.controller.geometry;  
  2.   
  3.   
  4. /** 
  5.  * @author BruceYang 
  6.  * 这个是对通用一次直线方程 A*x + B*y + C = 0 的封装~ 
  7.  * 本来封装的是斜截式,不过发现当斜率k不存在的时候会比较麻烦,因此该用一般式 
  8.  * 再个就是接着用一般式的演变方式 x + B/A*y + C/A = 0,但是考虑到可能存在x == 0 的情况,因此又舍弃~ 
  9.  *  
  10.  * 娘的,一般式还是他妈的无济于事啊,改回斜截式,多提供两个成员变量: 
  11.  * 一个boolean表示k是否存在,一个额外的float表示k不存在的时候直线方程 x=***, *** 等于多少~ 
  12.  */  
  13. public class CGLine {  
  14.     // 特别声明为public类型,免得到时候访问的时候麻烦,到时候直接点就行了  
  15.     private boolean kExists;    // 大部分情况下 k 都应该是存在的,因此提供一个 true 的默认值~  
  16.   
  17.     public double k = 77885.201314f;  
  18.     public double b = 13145.207788f;  
  19.     public double extraX = 52077.881314f;  
  20.       
  21.       
  22.     /** 
  23.      * 这是当 k 存在时的构造方法~ 
  24.      * @param k 
  25.      * @param b 
  26.      */  
  27.     public CGLine(double k, double b) {  
  28.         this.kExists = true;  
  29.         this.k = k;  
  30.         this.b = b;  
  31.     }  
  32.       
  33.     /** 
  34.      * 已知两点,求直线的方程~ 
  35.      * @param p1 
  36.      * @param p2 
  37.      */  
  38.     public CGLine(CGPoint p1, CGPoint p2) {  
  39.         if((p1.x - p2.x) != 0) {  
  40.             CGDbg.println("y = k*x + b, k exits!!");  
  41.             this.kExists = true;  
  42.             this.k = (p1.y - p2.y)/(p1.x - p2.x);  
  43.             this.b = (p1.y - p1.x * k);  
  44.         } else {  
  45.             CGDbg.println("y = k*x + b, k doesn't exists!!");  
  46.             // 如果走进这个分支,表示直线垂直于x轴,斜率不存在,保留k的默认值~  
  47.             this.kExists = false;  
  48.             this.extraX = p1.x;  
  49.         }  
  50.         CGDbg.print("过p1("+p1.x+", " +p1.y + "), p2("+p2.x+", "+p2.y+")两点的直线方程表达式为: ");  
  51.         if(kExists) {  
  52.             CGDbg.println("y = " + k + "*x + " + b);  
  53.         } else {  
  54.             CGDbg.println("x = " + extraX + "(垂直于x轴!)");  
  55.         }  
  56.     }  
  57.       
  58.     /** 
  59.      * 点斜式~ 
  60.      * @param p 某点 
  61.      * @param k 过该点的直线的斜率 
  62.      */  
  63.     public CGLine(double k, CGPoint p) {  
  64.         /** 
  65.          * (y-y') = k*(x-x') 
  66.          * 变形成斜截式为: 
  67.          * y = k*x + y' - k*x' 
  68.          * k = k, b = y'-k*x' 
  69.          */  
  70.         this.kExists = true;  
  71.         this.k = k;  
  72.         this.b = p.y - k * p.x;  
  73.     }  
  74.       
  75.     /** 
  76.      * 这是当 k 不存在时的构造方法~ 
  77.      * @param extraX 
  78.      */  
  79.     public CGLine(double extraX) {  
  80.         this.kExists = false;  
  81.         this.extraX = extraX;  
  82.     }  
  83.       
  84.     @Override  
  85.     public String toString() {  
  86.         return "Line.toString()方法被调用,y = k*x + b斜截式, k=" + this.k +   
  87.                 ", b=" + this.b +   
  88.                 ", kExists=" + this.kExists +   
  89.                 ", extraX=" + this.extraX;  
  90.     }  
  91.       
  92.     public boolean iskExists() {  
  93.         return kExists;  
  94.     }  
  95.     public void setkExists(boolean kExists) {  
  96.         this.kExists = kExists;  
  97.     }  
  98. }  

*****************************************************

[java] view plaincopy
  1. package org.bruce.vertices.controller.geometry;  
  2.   
  3. /** 
  4.  * @author Bruce Yang 
  5.  * 用于打印调试~ 
  6.  */  
  7. public class CGDbg {  
  8.     public static final boolean DEBUG_MODE = true;  
  9.       
  10.     // 方便进行调试信息的输出,开关~  
  11.     public static void println(Object info) {  
  12.         if(DEBUG_MODE) {              
  13.             System.out.println(info);  
  14.         }  
  15.     }  
  16.     public static void print(Object info) {  
  17.         if(DEBUG_MODE) {              
  18.             System.out.print(info);  
  19.         }  
  20.     }  
  21. }  

*****************************************************

[java] view plaincopy
  1. package org.bruce.vertices.controller.geometry;  
  2.   
  3. /** 
  4.  * @author BruceYang 
  5.  */  
  6. public class CGGeometryLib {  
  7.       
  8.     /** 
  9.      * @param p0    第一个点的坐标 
  10.      * @param p1    第二个点的坐标 
  11.      * @return      两个点之间的距离 
  12.      * 计算出两点之间的距离 
  13.      */  
  14.     public static double getDistanceBetween2Points(CGPoint p0, CGPoint p1) {  
  15.         double distance = Math.sqrt(Math.pow(p0.y - p1.y, 2) + Math.pow(p0.x - p1.x, 2));  
  16.         return distance;  
  17.     }  
  18.       
  19.     /** 
  20.      * @param p 
  21.      * @param l 
  22.      * @return      该方法用于获取某点在某条直线上的投影点的坐标 
  23.      */  
  24.     public static CGPoint getProjectivePoint(CGPoint p, CGLine l) {  
  25.         CGPoint target = null;  
  26.         if(l.iskExists()) {  
  27.             if(l.k != 0) {  
  28.                 CGLine newLine = new CGLine(-1/l.k, p.y -(-1/l.k)*p.x);  
  29.                 target = getCrossPoint(l, newLine);  
  30.             } else {    // 如果直线l的斜率存在且斜率的值为0,明显是一条平行于x轴的直线~  
  31.                 // 此时,点 p 到直线 l 的距离为:Math.abs(p.y-l.b)  
  32.                 target = new CGPoint(p.x, l.b);  
  33.             }  
  34.         } else {    // 如果直线l的斜率不存在,明显是一条垂直于x轴的直线~  
  35.             // 此时,点 p 到直线 l 的距离为:Math.abs(p.x-l.extraX)  
  36.             target = new CGPoint(l.extraX, p.y);  
  37.         }  
  38.         CGDbg.println("点 ("+p.x+", "+p.y+") 在直线:y="+l.k+"x+"+l.b+" 上的投影点为 ("+target.x+", "+target.y+")");  
  39.         return target;  
  40.     }  
  41.       
  42.     /** 
  43.      * 该方法用于求出两条直线的交点坐标 
  44.      * 这个方法是定制的,只有 l1, l2 均存在斜率 k 时方能使用(限制取消)~ 
  45.      * @param l1 
  46.      * @param l2 
  47.      * @return 
  48.      */  
  49.     public static CGPoint getCrossPoint(CGLine l1, CGLine l2) {  
  50. //      dbgPrintln("into the getCrossPoint, l1: " + l1);  
  51. //      dbgPrintln("into the getCrossPoint, l2: " + l2);  
  52.         double x, y;  
  53.         if(l1.iskExists() && l2.iskExists()) {  
  54.             x = (l2.b - l1.b) / (l1.k - l2.k);  
  55.             y = l1.k * x + l1.b;  
  56.         } else if(!l1.iskExists() && l2.iskExists()) {  
  57.             x = l1.extraX;  
  58.             y = l2.k * x + l2.b;  
  59.         } else if(l1.iskExists() && !l2.iskExists()) {  
  60.             x = l2.extraX;  
  61.             y = l1.k * x + l1.b;  
  62.         } else {  
  63.             // 两条直线的斜率都不存在?!,不可能发生的情况!!  
  64.             x = 0;  
  65.             y = 0;  
  66.         }  
  67.         CGDbg.println("getCrossPoint, CGPoint(x=" + x + ", y=" + y + ")");  
  68.         return new CGPoint(x, y);  
  69.     }  
  70.   
  71.     /** 
  72.      * @param args 
  73.      * 怎样判断是否符合要求?将过每组3个点中除开多边形顶点的两个点的直线方程求出来 
  74.      * 比较求出的 4 个 候选圆心点 哪个与此直线离的比较近,哪个就是符合要求的圆心点 
  75.      * 以下方法用于获取离特定直线距离最近的一个点(目前只支持斜率k存在的直线,以后慢慢扩充)! 
  76.      * 要得到距离特定直线距离最远的一个点只要稍作改动即可! 
  77.      */  
  78.     public static CGPoint getNearestPoint(CGPoint[] points, CGLine line) {  
  79.         double minDistance = 0;  
  80.         int minIndex = 0;  
  81.         if(line.iskExists()) {  
  82.             // 直线斜率存在的分支~  
  83.             for(int i = 0; i < points.length; ++ i) {  
  84.                 CGPoint p = points[i];  
  85.                 double d = Math.abs(line.k*p.x-p.y+line.b)/Math.sqrt(Math.pow(line.k,2)+1);  
  86.                 if(i == 0) {  
  87.                     // 赋予初值,不然 minDistance 的值就为 0 了~  
  88.                     minDistance = d;  
  89.                 }  
  90.                 if(d < minDistance) {  
  91.                     minDistance = d;  
  92.                     minIndex = i;  
  93.                 }  
  94.             }  
  95.         } else {  
  96.             // 直线斜率不存在的分支(亦即直线垂直于 x 轴)~  
  97.             for(int i = 0; i < points.length; ++ i) {  
  98.                 CGPoint p = points[i];  
  99.                 double d = Math.abs(p.x - line.extraX);  
  100.                 if(i == 0) {  
  101.                     // 赋予初值,不然minDistance的值就为0了~  
  102.                     minDistance = d;  
  103.                 }  
  104.                 if(d < minDistance) {  
  105.                     minDistance = d;  
  106.                     minIndex = i;  
  107.                 }  
  108.             }  
  109.         }  
  110.         CGPoint dest = points[minIndex];  
  111.         CGDbg.println("即将离开chooseOne()方法,圆心点为:("+dest.x+", "+dest.y+")");  
  112.         return dest;  
  113.     }  
  114.       
  115.     /** 
  116.      * 获取传入两点的中点~ 
  117.      * @param p1 
  118.      * @param p2 
  119.      * @return 
  120.      */  
  121.     public static CGPoint getMiddlePoint(CGPoint p1, CGPoint p2) {  
  122.         return new CGPoint((p1.x + p2.x) / 2.0f, (p1.y + p2.y) / 2.0f);  
  123.     }  
  124.       
  125.     /** 
  126.      * 封装一下 Math 的 pow 、sqrt 方法,调用起来方便一些~ 
  127.      * @param d1 
  128.      * @param d2 
  129.      * @return 
  130.      */  
  131.     public static double pow(double d1, double d2) {  
  132.         return Math.pow(d1, d2);  
  133.     }  
  134.     public static double sqrt(double d) {  
  135.         return Math.sqrt(d);  
  136.     }  
  137.     public static double sin(double theta) {  
  138.         return Math.sin(theta);  
  139.     }  
  140.     public static double cos(double theta) {  
  141.         return Math.cos(theta);  
  142.     }  
  143.       
  144.     /** 
  145.      * 传入线段的两个端点,获取中点,以该中点为圆心做半径为 radius 的圆, 
  146.      * 经过线段中点做线段的垂线,返回垂线与圆的两个交点~ 
  147.      * Objective-C 里面的结果有点儿问题,不知道是什么原因,来java 里面碰碰有运气~ 
  148.      * @param p1        线段端点1 
  149.      * @param p2        线段端点2 
  150.      * @param radius    圆半径 
  151.      * @return          线段中垂线与圆的两个交点~ 
  152.      */  
  153.     public static CGPoint[] getWhatIWanted(CGPoint p1, CGPoint p2, double radius) {  
  154.         CGPoint[] target = new CGPoint[2];  
  155.         CGPoint pMiddle = getMiddlePoint(p1, p2);  
  156. //      double segLength = getDistanceBetween2Points(p1, p2);  
  157.           
  158.         CGLine l1 = new CGLine(p1, p2);  
  159.         if(l1.iskExists()) {  
  160.             if(l1.k != 0) {  
  161.                 double kOfNewLine = -1 / l1.k;  
  162.                 CGLine newLine = new CGLine(kOfNewLine, pMiddle);  
  163.                   
  164.                 // 经过数学运算,得出二元一次方程组的表达式  
  165.                 double A = pow(newLine.k, 2) + 1;  
  166.                 double B = 2 * (newLine.k * newLine.b - newLine.k * pMiddle.y - pMiddle.x);  
  167.                 double C = pow(pMiddle.x, 2) + pow((newLine.b - pMiddle.y), 2) - pow(radius, 2);  
  168.                 double delta = pow(B, 2) - 4 * A * C;  
  169.                   
  170.                 if(delta < 0) {    // 经实践检验有一定几率走入该分支,必须做特殊化处理~  
  171.                     // 2012。04。28。20。01,精度不够所致,换成double后无该情况出现~  
  172.                     CGDbg.println("竟然会无解,他妈的怎么回事儿啊!");  
  173.                     target[0] = new CGPoint(pMiddle.x, pMiddle.y - radius);  
  174.                     target[1] = new CGPoint(pMiddle.x, pMiddle.y + radius);  
  175.                 } else {  
  176.                     double x1 = (-B + sqrt(delta)) / (2 * A);  
  177.                     double y1 = newLine.k * x1 + newLine.b;  
  178.                     target[0] = new CGPoint(x1, y1);  
  179.                       
  180.                     double x2 = (-B - sqrt(delta)) / (2 * A);  
  181.                     double y2 = newLine.k * x2 + newLine.b;  
  182.                     target[1] = new CGPoint(x2, y2);  
  183.                 }  
  184.             } else {  
  185.                 target[0] = new CGPoint(pMiddle.x, pMiddle.y - radius);  
  186.                 target[1] = new CGPoint(pMiddle.x, pMiddle.y + radius);  
  187.             }  
  188.         } else {  
  189.             target[0] = new CGPoint(pMiddle.x - radius, pMiddle.y);  
  190.             target[1] = new CGPoint(pMiddle.x + radius, pMiddle.y);  
  191.         }  
  192.         System.out.println("target[0] 距离中点的距离为:" + getDistanceBetween2Points(target[0], pMiddle));  
  193.         System.out.println("target[1] 距离中点的距离为:" + getDistanceBetween2Points(target[1], pMiddle));  
  194.         return target;  
  195.     }  
  196.       
  197.     /** 
  198.      * 测试实用性,测试结果如下: 
  199.      * 之前用 float 类型的时候,每隔 1 度测试一次,共测试一个圆周,无解的情况出现一次 
  200.      * 每隔 1 度测试一次, 共测试一个圆周,无解的情况无。 
  201.      * 每隔 0.5 度测试一次,共测试一个圆周,无解的情况只出现一次 
  202.      * @param args 
  203.      */  
  204.     public static void main(String[] args) {  
  205.         double currentRadian = 0;  
  206.         double deltaRadian = Math.PI / 360;  
  207.         double bigRadius = 50;  
  208.         double smallRadius = 20;  
  209.         CGPoint origin = new CGPoint(00); // 原点~  
  210.         CGPoint tail = null;    // tail 是尾巴、末梢的意思~  
  211.         for(int i = 0; i < 720; ++ i) {  
  212.             System.out.println(" -- 第 "+ i + "度!");  
  213.             tail = new CGPoint(bigRadius*cos(currentRadian), bigRadius*sin(currentRadian));  
  214.             currentRadian += deltaRadian;  
  215.             getWhatIWanted(origin, tail, smallRadius);  
  216.         }  
  217.     }  
  218. }  
原创粉丝点击