def visitArc(self, geometry): #self.logger.logMessageString("Arc feature", fmeobjects.FME_INFORM) # Circle through the three points is the problem of finding # the circumscribed circle about a triangle. # We use here cartesian eguations from page # https://en.wikipedia.org/wiki/Circumscribed_circle EPSILON_SQLMM = 0.00000001 startPoint = geometry.getStartPoint().getXYZ() endPoint = geometry.getEndPoint().getXYZ() threePointLine = geometry.getAsLineFixedArcSamples(3) midPoint = threePointLine.getPointAt(1).getXYZ() p1x=startPoint[0] p1y=startPoint[1] p1z=startPoint[2] p2x=midPoint[0] p2y=midPoint[1] p2z=midPoint[2] p3x=endPoint[0] p3y=endPoint[1] p3z=endPoint[2] #p1z = 100 #p2z = 50 #p3z = 200 dx21 = p2x - p1x dy21 = p2y - p1y dx31 = p3x - p1x dy31 = p3y - p1y # Special handling for circle # In PostGIS, if endpoints of the arc are same, then its geometry is # circle and midpoint means farthest opposite point of the circle. # We dont know the direction of the stroked line, but because circle # shouldn't never be part of the path, then it is not a problem. if ( abs(dx31) < EPSILON_SQLMM and abs(dy31) < EPSILON_SQLMM ): # If all the points are the same, then geometry is circle with zero radius # and we can return zero length line as its stroked geometry if ( abs(dx21) < EPSILON_SQLMM and abs(dy21) < EPSILON_SQLMM ): self.logger.logMessageString("Arc points were all the same. Returning zero length line as stroked geometry", fmeobjects.FME_INFORM) newGeometry = fmeobjects.FMELine([startPoint, endPoint]) if not geometry.is3D(): newGeometry.force2D() return newGeometry # d < 0 means clockwise sweeping angle # d > 0 means counterclockwise sweeping angle # We use here for circle counterclockwise sweeping angle d = 1 # Circle origo X coordinate cx = p1x + (p2x - p1x) / 2.0 # Circle origo Y coordinate cy = p1y + (p2y - p1y) / 2.0; # Circle radius r = math.sqrt((p1x - cx)**2 + (p1y - cy)**2) # Angle counterclockwise from x-axis to p1 a1 = math.atan2(p1y - cy, p1x - cx) # Angle counterclockwise from x-axis to p2 a2 = a1 + math.copysign(1, d) * math.pi # Angle counterclockwise from x-axis to p3 a3 = a1 + math.copysign(1, d) * 2 * math.pi # Geometry is not a circle but an arc, so we are looking for radius and centroid of the arc else: h21 = dx21**2 + dy21**2 h31 = dx31**2 + dy31**2 # 2 * |Cross product| # d < 0 means clockwise sweeping angle # d > 0 means counterclockwise sweeping angle d = 2 * (dx21 * dy31 - dx31 * dy21) # |Cross product| can be zero only when all points are in same line. # Replacing arc with straight line should be correct stroked geometry here. # To be strict, if midpoint is same as startpoint or endpoint, then arc is not defined at all. # That kind of geometries shouldn't exist, but here we return straight line for them also. if abs(d) < EPSILON_SQLMM: self.logger.logMessageString("Arc points were in same line. Returning straight line as stroked geometry", fmeobjects.FME_INFORM) newGeometry = fmeobjects.FMELine([startPoint, endPoint]) if not geometry.is3D(): newGeometry.force2D() return newGeometry # Arc origo X coordinate cx = p1x + (h21 * dy31 - h31 * dy21) / d # Arc origo Y coordinate cy = p1y - (h21 * dx31 - h31 * dx21) / d # Arc radius r = math.sqrt((p1x - cx)**2 + (p1y - cy)**2) # Angle counterclockwise from x-axis to p1 a1 = math.atan2(p1y - cy, p1x - cx) # Angle counterclockwise from x-axis to p2 a2 = math.atan2(p2y - cy, p2x - cx) # Make needed correction on a2. Sign of d gives sweep direction. a2 = a2 + (math.copysign(1, d) - math.copysign(1, a2 - a1)) * math.pi # Angle counterclockwise from x-axis to p3 a3 = math.atan2(p3y - cy, p3x - cx) # Make needed correction on a3. Sign of d gives sweep direction. a3 = a3 + (math.copysign(1, d) - math.copysign(1, a3 - a1)) * math.pi # Radius, centroid and angles are now calculated to geometry. # Then we start stroking the geometry angleArc = a3 - a1 toleranceVertices = 0 verticesCount = 0 StrokingTolerance = 0.0 # If user gave stroking tolerance, we calculate maximum angle of one stroke sector. # from eguation: cos(angle / 2) = (radius - tolerance) / radius # We calculate also minimum count of vertices for whole arc that satisfies tolerance. if self.__arcStrokingTolerance > 0.0: StrokingTolerance = self.__arcStrokingTolerance if StrokingTolerance > 2 * r: StrokingTolerance = 2 * r angleSectorFromTolerance = 2 * math.acos((r - StrokingTolerance) / r) toleranceVertices = int(abs(angleArc / angleSectorFromTolerance)) # If user gave exact amout of vertices to use, # we use it or we use vertice count calculated from tolerance, # wich ever is bigger so meaning more accurate if self.__arcStrokingVertices > toleranceVertices: verticesCount = self.__arcStrokingVertices else: verticesCount = toleranceVertices # Now calculate exact angle of one stroke sector from count of vertices angleSector = angleArc / (verticesCount + 1) calculateZ = True if (p1z == p2z) and (p1z == p3z): calculateZ = False newPointZ = p1z # Calculate vertice points and create line from them pointList = [(p1x, p1y, p1z)] for i in range(1, verticesCount + 1): currentAngle = a1 + i * angleSector newPointX = cx + r * math.cos(currentAngle) newPointY = cy + r * math.sin(currentAngle) if calculateZ: if math.copysign(1, a2 - currentAngle) == math.copysign(1, d): newPointZ = p1z + (p2z - p1z) * (currentAngle - a1) / (a2 - a1) else: newPointZ = p2z + (p3z - p2z) * (currentAngle - a2) / (a3 - a2) pointList.append((newPointX, newPointY, newPointZ)) pointList.append((p3x, p3y, p3z)) newGeometry = fmeobjects.FMELine(pointList) if not geometry.is3D(): newGeometry.force2D() return newGeometry