Bulge

Ура, я почти закончил! На свободу с чистой совестью!

Не так это и сложно, спустя 20 лет вспомнить тригонометрию.

Есть в чертежах (DWG) AutoCAD такой примитив как polyline, иногда LWPolyline. Ничего особенного, набор точек, соединенных отрезками. Подвох в том, что у любого отрезка (определяемого двумя соседними точками) может быть ненулевым такой параметр как bulge. Эту хрень придумали чтоб мне было чем заняться чтобы можно было в полилинию встраивать не только прямые отрезки, но и дуги (arc).

Вот как описывает bulge один почтипоэт:

Bulges are something that women have (mostly to please the opposite sex it seems) and something that guys try to get by placing socks in strategic places. At least until they get older. Which is the time they tend to develop bulges in not so strategic places. In other words: bulges are all about curvature.

...

In AutoCAD's online help reference, it says about bulges for polylines:

The bulge is the tangent of 1/4 of the included angle for the arc between the selected vertex and the next vertex in the polyline's vertex list. A negative bulge value indicates that the arc goes clockwise from the selected vertex to the next vertex. A bulge of 0 indicates a straight segment, and a bulge of 1 is a semicircle.

afralisp.net/archive/lisp/Bulges

Предупреждаю, если полезете по сцылке, рекомендую отключить в браузере использование стилей. В моем ФФ4 видно белый текст на белом фоне, если со стилями.

А мне надо эти полилинии с булжами преобразовать в обычные полилинии, где кроме прямых отрезков нет ничего. Соответственно, дуги надо аппроксимировать отрезками. Правда, сначала надо превратить две точки и число (bulge) в геометрическую фигуру — дугу. Но оказалось не слишком все тяжко. Полдня поисков, полдня воспоминаний школьного курса тригонометрии и полдня передирания формул из ёкселя в Python.

Черновой вариант превращения сегмента полилинии с булжем в обычную полилинию готов. 

UPD

Не готов, ибо начиная с центра дуги, считает неправильно. Когда будет правильно - обновлю пост. А пока так. 

UPD2

Вот теперь правильно. Проверено отладкой на сегментах в разных квадрантах и с разным направлением.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# (c) Valik mailto:vasnake@gmail.com

import sys, os, math

def autoLispAngle(x1,y1, x2,y2):
 return math.atan2(y2-y1, x2-x1)

def polar(x1,y1, phi, dist):
 x = x1 + dist * math.cos(phi)
 y = y1 + dist * math.sin(phi)
 return (x,y)

def sign(n):
 if n < 0.0: return -1.0
 return 1.0

def testBulge(x1, y1, x2, y2, bulge, sublen=0.01):
 ''' AutoCAD polyline arc segment approximation.
 Bulge to arc to line segments (facets)
 http://www.cadtutor.net/forum/showthread.php?51511-Points-along-a-lwpoly-arc&
 att\FacetBulge_rev1.zip\FacetBulge\BulgeCalc_Rev1.xls - formulas with errors, which eliminate with:
 http://www.afralisp.net/archive/lisp/Bulges1.htm
 '''
 print 'calcBulge: p1 [%s], p2 [%s], bulge [%s], sublen [%s]' % ((x1,y1), (x2,y2), bulge, sublen)
 res = {}
# midpoint (F9, G9) # arc midpoint [(12.41925, 13.5352)]; // 16.73, 11.77
 mx1 = (x1 + x2) / 2.0
 my1 = (y1 + y2) / 2.0
 s1 = 'chord midpoint [%s]' % ((mx1, my1),)
 res['chordmidpoint'] = (mx1, my1)
# angle (K13) # angle/2 (K14)
 angle = math.atan(bulge) * 4.0
 angleDeg = angle * (180.0 / math.pi)
 s2 = 'arc angle [%0.5f] rad, [%0.5f] deg' % (angle, angleDeg)
 res['angleRad'] = angle
 res['angleDeg'] = angleDeg
# dist (F5)
 dist = math.sqrt((x2-x1)**2 + (y2-y1)**2)
 s3 = 'chord len [%0.5f]' % dist
 res['chordLen'] = dist
# sagitta length (http://en.wikipedia.org/wiki/Sagitta_%28geometry%29)
 sagitta = dist/2.0 * bulge
# radius (K15) # r = ((dist/2.0)**2+sagitta**2)/2.0*sagitta
 radius = (dist/2.0) / math.sin( abs(angle/2.0) )
 if radius == 0.0: radius = 0.000000000000001
 s4 = 'radius [%0.5f]' % radius
 res['radius'] = radius
# arc length (F14) l = 2*pi*r
 alen = abs(radius * angle)
 s5 = 'arc length [%0.5f]' % alen
 res['arcLen'] = alen
# center (F19, G19)
 t = '''
wrong algo:
 k5 = ( (math.sqrt(radius**2 - (dist / 2.0)**2))*2 ) / dist
 cx = mx1 + (((y1-y2)/2.0) * k5 * sign(bulge))
 cy = my1 + (((x2-x1)/2.0) * k5 * sign(bulge))
 s6 = 'arc center [%s]' % ((cx, cy),)
 res['center'] = (cx, cy)
good algo:
 (setq bulge 2.5613 p1 (list 11.7326 11.8487) p2 (list 13.1059 15.2217) r 2.68744
   theta (* 4.0 (atan (abs bulge)))
   gamma (/ (- pi theta) 2.0)
   phi (+ (angle p1 p2) gamma)
   p (polar p1 phi r)
 )
 '''
 theta = 4.0 * math.atan(abs(bulge))
 gamma = (math.pi - theta) / 2.0
 phi = autoLispAngle(x1,y1, x2,y2) + gamma * sign(bulge)
 cx,cy = polar(x1,y1, phi, radius)
 s6 = 'arc center [%s]' % ((cx, cy),) # 14.2498, 12.7899
 res['center'] = (cx, cy)
# start, end angle (G21, G22)
 startAngle = math.acos((x1 - cx) / radius)
 if not sign(y1 - cy) > 0:
  startAngle = (2.0 * math.pi) - startAngle
 endAngle = startAngle + angle
 s7 = 'start, end arc angles [%s]' % ((startAngle, endAngle),)
 res['seAngles'] = (startAngle, endAngle)
# subangle (F27), numsub (# of Divisions K26)
 if sublen <= 0.0: sublen = alen / 10.0
 numsub = round(alen/sublen, 0)
 if numsub < 2:
  numsub = 2.0
 subangle = angle / numsub
 s8 = 'numsub, subangle [%s]' % ((numsub, subangle),)
# length of subarc L27
 realSublen = abs(2 * (math.sin(subangle/2.0) * radius) )
 s9 = 'real sublen [%0.5f]' % realSublen
# sub points
 currangle = startAngle + (subangle/2.0) + (math.pi/2.0 * sign(bulge))
 sx = x1
 sy = y1
 listPoints = [(x1,y1)]
 for cnt in range(int(numsub-1)):
  if not cnt == 0:
   currangle = currangle + subangle
  sx = sx + (realSublen * math.cos(currangle))
  sy = sy + (realSublen * math.sin(currangle))
  listPoints.append((sx, sy))
 listPoints.append((x2, y2))
 res['points'] = listPoints
 print 'doBulge done [%s; %s; %s; %s; %s; %s; %s; %s; %s; subpoints %s]' % (s1,s2,s3,s4,s5,s6,s7,s8,s9,listPoints)
 return res
#def testBulge(self, x1, y1, x2, y2, bulge, sublen=10):

def main():
 testBulge(13.0, 9.0, 21.68, 33.65, -0.45, 7)
 #~ testBulge(11.7326, 11.8487, 13.1059, 15.2217, 2.5613, 1)
 #~ testBulge(5.9228, 12.0274, 8.1062, 22.8887, 2.2893, 3)
 #~ testBulge(24.2884, 10.3276, 27.3493, 14.9170, -4.5373, 3)
 #~ testBulge(-34.8952, -21.6100, -32.7117, -10.7486, 2.2893, 5)
 #~ testBulge(-16.5296, -23.3098, -13.4687, -18.7203, -4.5373, 5)
 #~ testBulge(10.7004, -22.2329, 12.8839, -11.3716, 2.2893, 2)
 #~ testBulge(29.0660, -23.9327, 32.1270, -19.3433, -4.5373, 2)
 #~ testBulge(-34.2720, 12.0274, -32.0886, 22.8887, 2.2893, 2)
 #~ testBulge(-15.9064, 10.3276, -12.8455, 14.9170, -4.5373, 2)
 #~ testBulge(8.0555, 5.0696, -6.8401, -6.4998, -0.3324, 2)
 #~ testBulge(-3.6195, 4.3654, 6.0425, -6.4998, 2.1734, 2)
 #~ testBulge(6.0425, -6.4998, -3.6195, 4.3654, -2.1734, 2)
 #~ testBulge(1.0, 0.0, 0.0, 1.0, 1, 0.5)
 #~ testBulge(0, 1, -1, 0, 1, 0.5)
 #~ testBulge(1.0, 0.0, 0.99, -0.044, 44.53, 0.3)
 #~ testBulge(0.0, 0.0, 0.0, 0.0, 45, 0.3)

import time, traceback
print time.strftime('%Y-%m-%d %H:%M:%S')
if __name__ == "__main__":
 try:
  main()
 except Exception, e:
  if type(e).__name__ == 'COMError': print 'COM Error, msg [%s]' % e
  else:
   print 'Error, program failed:'
   traceback.print_exc(file=sys.stderr)
print time.strftime('%Y-%m-%d %H:%M:%S')

Надо будет подчистить, убрать накопление ошибки при аппроксимации большим количеством сегментов, обобщить, чтобы и просто дуги с окружностями (circle) можно было в полилинии превращать. Ну и отладить эту байду. Еще полдня — день.

Спасибо Autodesk-у за то, что у меня есть работа, причем достаточно интересная.

Это была реклама Autodesk и его софтины AutoCAD :)