diff --git a/src/Mod/Draft/CMakeLists.txt b/src/Mod/Draft/CMakeLists.txt index f083e45631..69c2d5263d 100644 --- a/src/Mod/Draft/CMakeLists.txt +++ b/src/Mod/Draft/CMakeLists.txt @@ -43,6 +43,7 @@ SET (Draft_geoutils draftgeoutils/linear_algebra.py draftgeoutils/cuboids.py draftgeoutils/circles.py + draftgeoutils/circles_apollonius.py ) SET(Draft_tests diff --git a/src/Mod/Draft/DraftGeomUtils.py b/src/Mod/Draft/DraftGeomUtils.py index 55a8a71279..9e986d3459 100644 --- a/src/Mod/Draft/DraftGeomUtils.py +++ b/src/Mod/Draft/DraftGeomUtils.py @@ -385,167 +385,13 @@ from draftgeoutils.arcs import arcFrom2Pts #############################33 to include +from draftgeoutils.circles_apollonius import outerSoddyCircle +from draftgeoutils.circles_apollonius import innerSoddyCircle - - -def outerSoddyCircle(circle1, circle2, circle3): - """Compute the outer soddy circle for three tightly packed circles.""" - if (geomType(circle1) == "Circle") and (geomType(circle2) == "Circle") \ - and (geomType(circle3) == "Circle"): - # Original Java code Copyright (rc) 2008 Werner Randelshofer - # Converted to python by Martin Buerbaum 2009 - # http://www.randelshofer.ch/treeviz/ - # Either Creative Commons Attribution 3.0, the MIT license, or the GNU Lesser General License LGPL. - - A = circle1.Curve.Center - B = circle2.Curve.Center - C = circle3.Curve.Center - - ra = circle1.Curve.Radius - rb = circle2.Curve.Radius - rc = circle3.Curve.Radius - - # Solution using Descartes' theorem, as described here: - # http://en.wikipedia.org/wiki/Descartes%27_theorem - k1 = 1 / ra - k2 = 1 / rb - k3 = 1 / rc - k4 = abs(k1 + k2 + k3 - 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1)) - - q1 = (k1 + 0j) * (A.x + A.y * 1j) - q2 = (k2 + 0j) * (B.x + B.y * 1j) - q3 = (k3 + 0j) * (C.x + C.y * 1j) - - temp = ((q1 * q2) + (q2 * q3) + (q3 * q1)) - q4 = q1 + q2 + q3 - ((2 + 0j) * cmath.sqrt(temp) ) - - z = q4 / (k4 + 0j) - - # If the formula is not solvable, we return no circle. - if (not z or not (1 / k4)): - return None - - X = -z.real - Y = -z.imag - print("Outer Soddy circle: " + str(X) + " " + str(Y) + "\n") # Debug - - # The Radius of the outer soddy circle can also be calculated with the following formula: - # radiusOuter = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 - 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3)))) - circ = Part.Circle(Vector(X, Y, A.z), norm, 1 / k4) - return circ - - else: - print("debug: outerSoddyCircle bad parameters!\n") - # FreeCAD.Console.PrintMessage("debug: outerSoddyCircle bad parameters!\n") - return None - - -def innerSoddyCircle(circle1, circle2, circle3): - """Compute the inner soddy circle for three tightly packed circles.""" - if (geomType(circle1) == "Circle") and (geomType(circle2) == "Circle") \ - and (geomType(circle3) == "Circle"): - # Original Java code Copyright (rc) 2008 Werner Randelshofer - # Converted to python by Martin Buerbaum 2009 - # http://www.randelshofer.ch/treeviz/ - - A = circle1.Curve.Center - B = circle2.Curve.Center - C = circle3.Curve.Center - - ra = circle1.Curve.Radius - rb = circle2.Curve.Radius - rc = circle3.Curve.Radius - - # Solution using Descartes' theorem, as described here: - # http://en.wikipedia.org/wiki/Descartes%27_theorem - k1 = 1 / ra - k2 = 1 / rb - k3 = 1 / rc - k4 = abs(k1 + k2 + k3 + 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1)) - - q1 = (k1 + 0j) * (A.x + A.y * 1j) - q2 = (k2 + 0j) * (B.x + B.y * 1j) - q3 = (k3 + 0j) * (C.x + C.y * 1j) - - temp = ((q1 * q2) + (q2 * q3) + (q3 * q1)) - q4 = q1 + q2 + q3 + ((2 + 0j) * cmath.sqrt(temp) ) - - z = q4 / (k4 + 0j) - - # If the formula is not solvable, we return no circle. - if (not z or not (1 / k4)): - return None - - X = z.real - Y = z.imag - print("Outer Soddy circle: " + str(X) + " " + str(Y) + "\n") # Debug - - # The Radius of the inner soddy circle can also be calculated with the following formula: - # radiusInner = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 + 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3)))) - circ = Part.Circle(Vector(X, Y, A.z), norm, 1 / k4) - return circ - - else: - print("debug: innerSoddyCircle bad parameters!\n") - # FreeCAD.Console.PrintMessage("debug: innerSoddyCircle bad parameters!\n") - return None - - -def circleFrom3CircleTangents(circle1, circle2, circle3): - """ - http://en.wikipedia.org/wiki/Problem_of_Apollonius#Inversive_methods - http://mathworld.wolfram.com/ApolloniusCircle.html - http://mathworld.wolfram.com/ApolloniusProblem.html - """ - - if (geomType(circle1) == "Circle") and (geomType(circle2) == "Circle") \ - and (geomType(circle3) == "Circle"): - int12 = findIntersection(circle1, circle2, True, True) - int23 = findIntersection(circle2, circle3, True, True) - int31 = findIntersection(circle3, circle1, True, True) - - if int12 and int23 and int31: - if len(int12) == 1 and len(int23) == 1 and len(int31) == 1: - # Only one intersection with each circle. - # => "Soddy Circle" - 2 solutions. - # http://en.wikipedia.org/wiki/Problem_of_Apollonius#Mutually_tangent_given_circles:_Soddy.27s_circles_and_Descartes.27_theorem - # http://mathworld.wolfram.com/SoddyCircles.html - # http://mathworld.wolfram.com/InnerSoddyCenter.html - # http://mathworld.wolfram.com/OuterSoddyCenter.html - - r1 = circle1.Curve.Radius - r2 = circle2.Curve.Radius - r3 = circle3.Curve.Radius - outerSoddy = outerSoddyCircle(circle1, circle2, circle3) -# print(str(outerSoddy) + "\n") # Debug - - innerSoddy = innerSoddyCircle(circle1, circle2, circle3) -# print(str(innerSoddy) + "\n") # Debug - - circles = [] - if outerSoddy: - circles.append(outerSoddy) - if innerSoddy: - circles.append(innerSoddy) - return circles - - # @todo Calc all 6 homothetic centers. - # @todo Create 3 lines from the inner and 4 from the outer h. center. - # @todo Calc. the 4 inversion poles of these lines for each circle. - # @todo Calc. the radical center of the 3 circles. - # @todo Calc. the intersection points (max. 8) of 4 lines (through each inversion pole and the radical center) with the circle. - # This gives us all the tangent points. - else: - # Some circles are inside each other or an error has occurred. - return None - - else: - print("debug: circleFrom3CircleTangents bad parameters!\n") - # FreeCAD.Console.PrintMessage("debug: circleFrom3CircleTangents bad parameters!\n") - return None +from draftgeoutils.circles_apollonius import circleFrom3CircleTangents from draftgeoutils.linear_algebra import linearFromPoints diff --git a/src/Mod/Draft/draftgeoutils/circles_apollonius.py b/src/Mod/Draft/draftgeoutils/circles_apollonius.py new file mode 100644 index 0000000000..c52fbafcb1 --- /dev/null +++ b/src/Mod/Draft/draftgeoutils/circles_apollonius.py @@ -0,0 +1,222 @@ +# *************************************************************************** +# * Copyright (c) 2009, 2010 Yorik van Havre * +# * Copyright (c) 2009, 2010 Ken Cline * +# * * +# * This file is part of the FreeCAD CAx development system. * +# * * +# * This program is free software; you can redistribute it and/or modify * +# * it under the terms of the GNU Lesser General Public License (LGPL) * +# * as published by the Free Software Foundation; either version 2 of * +# * the License, or (at your option) any later version. * +# * for detail see the LICENCE text file. * +# * * +# * FreeCAD is distributed in the hope that it will be useful, * +# * but WITHOUT ANY WARRANTY; without even the implied warranty of * +# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +# * GNU Library General Public License for more details. * +# * * +# * You should have received a copy of the GNU Library General Public * +# * License along with FreeCAD; if not, write to the Free Software * +# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * +# * USA * +# * * +# *************************************************************************** +"""Provides various functions for Appollonius and Soddy circle operations. + +The 'problem of Appollonius' consists of finding a circle that is +simultaneously tangent to three circles on a plane. + +- http://en.wikipedia.org/wiki/Problem_of_Apollonius#Inversive_methods +- http://mathworld.wolfram.com/ApolloniusCircle.html +- http://mathworld.wolfram.com/ApolloniusProblem.html + +If each circle only has one intersection with each other, the problem is that +of finding a 'Soddy circle', which has two solutions, and inner circle +and an outer circle. + +- http://en.wikipedia.org/wiki/Problem_of_Apollonius#Mutually_tangent_given_circles:_Soddy.27s_circles_and_Descartes.27_theorem +- http://mathworld.wolfram.com/SoddyCircles.html +- http://mathworld.wolfram.com/InnerSoddyCenter.html +- http://mathworld.wolfram.com/OuterSoddyCenter.html +""" +## @package circles_apollonius +# \ingroup DRAFTGEOUTILS +# \brief Provides various functions for Appollonius and Soddy circles. + +import cmath +import math +import lazy_loader.lazy_loader as lz + +import FreeCAD as App +import DraftVecUtils + +from draftgeoutils.general import geomType, vec, NORM +from draftgeoutils.intersections import findIntersection + + +# Delay import of module until first use because it is heavy +Part = lz.LazyLoader("Part", globals(), "Part") + + +def outerSoddyCircle(circle1, circle2, circle3): + """Compute the outer soddy circle for three tightly packed circles. + + Original Java code Copyright (rc) 2008 Werner Randelshofer + Converted to python by Martin Buerbaum 2009 + http://www.randelshofer.ch/treeviz/ + Either Creative Commons Attribution 3.0, the MIT license, + or the GNU Lesser General License LGPL. + """ + if (geomType(circle1) != "Circle" + or geomType(circle2) != "Circle" + or geomType(circle3) != "Circle"): + print("debug: outerSoddyCircle bad parameters! Must be circles.") + return None + + A = circle1.Curve.Center + B = circle2.Curve.Center + C = circle3.Curve.Center + + ra = circle1.Curve.Radius + rb = circle2.Curve.Radius + rc = circle3.Curve.Radius + + # Solution using Descartes' theorem, as described here: + # http://en.wikipedia.org/wiki/Descartes%27_theorem + k1 = 1 / ra + k2 = 1 / rb + k3 = 1 / rc + k4 = abs(k1 + k2 + k3 - 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1)) + + q1 = (k1 + 0j) * (A.x + A.y * 1j) + q2 = (k2 + 0j) * (B.x + B.y * 1j) + q3 = (k3 + 0j) * (C.x + C.y * 1j) + + temp = ((q1 * q2) + (q2 * q3) + (q3 * q1)) + q4 = q1 + q2 + q3 - ((2 + 0j) * cmath.sqrt(temp)) + + z = q4 / (k4 + 0j) + + # If the formula is not solvable, we return no circle. + if not z or not (1 / k4): + return None + + X = -z.real + Y = -z.imag + print("Outer Soddy circle: " + str(X) + " " + str(Y) + "\n") # Debug + + # The Radius of the outer soddy circle can also be calculated + # with the following formula: + # radiusOuter = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 - 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3)))) + circ = Part.Circle(App.Vector(X, Y, A.z), NORM, 1 / k4) + + return circ + + +def innerSoddyCircle(circle1, circle2, circle3): + """Compute the inner soddy circle for three tightly packed circles. + + Original Java code Copyright (rc) 2008 Werner Randelshofer + Converted to python by Martin Buerbaum 2009 + http://www.randelshofer.ch/treeviz/ + """ + if (geomType(circle1) != "Circle" + and geomType(circle2) != "Circle" + and geomType(circle3) != "Circle"): + print("debug: innerSoddyCircle bad parameters! Must be circles.") + return None + + A = circle1.Curve.Center + B = circle2.Curve.Center + C = circle3.Curve.Center + + ra = circle1.Curve.Radius + rb = circle2.Curve.Radius + rc = circle3.Curve.Radius + + # Solution using Descartes' theorem, as described here: + # http://en.wikipedia.org/wiki/Descartes%27_theorem + k1 = 1 / ra + k2 = 1 / rb + k3 = 1 / rc + k4 = abs(k1 + k2 + k3 + 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1)) + + q1 = (k1 + 0j) * (A.x + A.y * 1j) + q2 = (k2 + 0j) * (B.x + B.y * 1j) + q3 = (k3 + 0j) * (C.x + C.y * 1j) + + temp = ((q1 * q2) + (q2 * q3) + (q3 * q1)) + q4 = q1 + q2 + q3 + ((2 + 0j) * cmath.sqrt(temp)) + + z = q4 / (k4 + 0j) + + # If the formula is not solvable, we return no circle. + if (not z or not (1 / k4)): + return None + + X = z.real + Y = z.imag + print("Outer Soddy circle: " + str(X) + " " + str(Y) + "\n") # Debug + + # The Radius of the inner soddy circle can also be calculated + # with the following formula: + # radiusInner = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 + 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3)))) + circ = Part.Circle(App.Vector(X, Y, A.z), NORM, 1 / k4) + + return circ + + +def circleFrom3CircleTangents(circle1, circle2, circle3): + """Return the circle that is tangent to three other circles. + + This problem is called the 'Problem of Appollonius'. + + A special case is that of 'Soddy circles'. + + To Do + ----- + Currently not all possible solutions are found, only the Soddy circles. + + * Calc all 6 homothetic centers. + * Create 3 lines from the inner and 4 from the outer h. center. + * Calc. the 4 inversion poles of these lines for each circle. + * Calc. the radical center of the 3 circles. + * Calc. the intersection points (max. 8) of 4 lines (through each + inversion pole and the radical center) with the circle. + * This gives us all the tangent points. + """ + if (geomType(circle1) != "Circle" + and geomType(circle2) != "Circle" + and geomType(circle3) == "Circle"): + print("debug: circleFrom3CircleTangents bad input! Must be circles") + return None + + int12 = findIntersection(circle1, circle2, True, True) + int23 = findIntersection(circle2, circle3, True, True) + int31 = findIntersection(circle3, circle1, True, True) + + if int12 and int23 and int31: + if len(int12) == 1 and len(int23) == 1 and len(int31) == 1: + # If only one intersection with each circle, Soddy circle + + # r1 = circle1.Curve.Radius + # r2 = circle2.Curve.Radius + # r3 = circle3.Curve.Radius + outerSoddy = outerSoddyCircle(circle1, circle2, circle3) + # print(str(outerSoddy)) # Debug + + innerSoddy = innerSoddyCircle(circle1, circle2, circle3) + # print(str(innerSoddy)) # Debug + + circles = [] + if outerSoddy: + circles.append(outerSoddy) + if innerSoddy: + circles.append(innerSoddy) + return circles + + # Here the rest of the circles should be calculated + # ... + else: + # Some circles are inside each other or an error has occurred. + return None