UPDATE: My code here has a bug. Mr. Rafsanjani corrected it in the comments. Please refer to them below. Apologies - this blog is for beginners and for learning, so I've left the original code and the comment as is. CBT 21NOV2014
I had previously done a series of posts on polygon offset intended as a practical guide for accomplishing the task quickly. Some kind folks pointed out that I was making things considerably harder than necessary by using trigonometric functions when vector math would be easier and less error prone.
A coworker lent me his Java code that does polygon offset. I translated it into python (2.5) using the pyeuclid module:
import euclid as eu
import copy
OFFSET = 0.15
# coordinates
# PT 1
MONASTERY = [(1.1, 0.75),
# PT 2
(1.2, 1.95),
. . .
# PT 21
(1.1, 0.75)]
def scaleadd(origin, offset, vectorx):
"""
From a vector representing the origin,
a scalar offset, and a vector, returns
a Vector3 object representing a point
offset from the origin.
(Multiply vectorx by offset and add to origin.)
"""
multx = vectorx * offset
return multx + origin
def getinsetpoint(pt1, pt2, pt3):
"""
Given three points that form a corner (pt1, pt2, pt3),
returns a point offset distance OFFSET to the right
of the path formed by pt1-pt2-pt3.
pt1, pt2, and pt3 are two tuples.
Returns a Vector3 object.
"""
origin = eu.Vector3(pt2[0], pt2[1], 0.0)
v1 = eu.Vector3(pt1[0] - pt2[0],
pt1[1] - pt2[1], 0.0)
v1.normalize()
v2 = eu.Vector3(pt3[0] - pt2[0],
pt3[1] - pt2[1], 0.0)
v2.normalize()
v3 = copy.copy(v1)
v1 = v1.cross(v2)
v3 += v2
if v1.z < 0.0:
retval = scaleadd(origin, -OFFSET, v3)
else:
retval = scaleadd(origin, OFFSET, v3)
return retval
polyinset = []
lenpolygon = len(MONASTERY)
i = 0
poly = MONASTERY
while i < lenpolygon - 2:
polyinset.append(getinsetpoint(poly[i],
poly[i + 1], poly[i + 2]))
i += 1
polyinset.append(getinsetpoint(poly[-2],
poly[0], poly[1]))
polyinset.append(getinsetpoint(poly[0],
poly[1], poly[2]))
The result:
Hey thanks heaps for writing this up. You should have seen the crazy quadrant atan2 path I was going down!
ReplyDeleteCould you link in your post the Java code your coworker passed you on?
ReplyDeleteMiguel,
ReplyDeleteI'm sorry. I no longer have access to the Java code.
Carl T.
It's a pity. Thanks anyway.
ReplyDeleteDear Carl,
ReplyDeleteMany thanks for this nice code.
I think there is a small bug:
In "getinsetpoint", the vector v3 should be normalized before passing to "scaleadd".
Furthermore, the final offset is not as the prescribed OFFSET and the angle between vectors should be taken into account.
A possible solution could be:
def getinsetpoint(pt1, pt2, pt3,offset):
"""
Given three points that form a corner (pt1, pt2, pt3),
returns a point offset distance OFFSET to the right
of the path formed by pt1-pt2-pt3.
pt1, pt2, and pt3 are two tuples.
Returns a Vector3 object.
"""
origin = eu.Vector3(pt2[0], pt2[1], 0.0)
v1 = eu.Vector3(pt1[0] - pt2[0],
pt1[1] - pt2[1], 0.0)
v1.normalize()
v2 = eu.Vector3(pt3[0] - pt2[0],
pt3[1] - pt2[1], 0.0)
v2.normalize()
v3 = copy.copy(v1)
v1 = v1.cross(v2)
v3 += v2
v3.normalize()
cs=v3.dot(v2)
a1=cs*v2
a2=v3-a1
if cs>0:
alpha=sqrt(a2.magnitude_squared())
else:
alpha=-sqrt(a2.magnitude_squared())
if v1.z < 0.0:
retval = scaleadd(origin, -offset/alpha, v3)
else:
retval = scaleadd(origin, offset/alpha, v3)
return retval
@Ahmad Thank you for stopping by and for reviewing the code. I don't have gnuplot handy at the moment to see what would happen if I instituted the changes you suggest. I think what I did originally is just inspect the plot and felt that it looked OK. What I need to do is measure the actual offset in the code and report it. I probably won't get to this immediately, but I intend to take a second look at this.
ReplyDeleteAgain, thank you for the rigor with which you reviewed this code. Carl T.
Mr. Rafsanjani,
ReplyDeleteI finally got around to updating this with some data that illustrates the problem you pointed out:
http://pyright.blogspot.com/2014/11/polygon-offset-with-pyeuclid-revisited.html
Thanks so much for your help. Your attention to detail and willingness to point me in the right direction are appreciated. CBT