Skip to content
Permalink
dev
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
import math
import numpy as np
def closestDistanceBetweenLines(a0, a1, b0, b1, clampAll=False, clampA0=False, clampA1=False, clampB0=False, clampB1=False):
''' Given two lines defined by numpy.array pairs (a0,a1,b0,b1)
Return distance, the two closest points, and their average
'''
# If clampAll=True, set all clamps to True
if clampAll:
clampA0 = True
clampA1 = True
clampB0 = True
clampB1 = True
# Calculate denomitator
A = a1 - a0
B = b1 - b0
normA = np.linalg.norm(B)
if normA == 0:
print("HERE", b0, b1)
_A = A / np.linalg.norm(A)
_B = B / np.linalg.norm(B)
cross = np.cross(_A, _B);
denom = np.linalg.norm(cross) ** 2
# If denominator is 0, lines are parallel: Calculate distance with a projection
# and evaluate clamp edge cases
if (denom == 0):
d0 = np.dot(_A, (b0 - a0))
d = np.linalg.norm(((d0 * _A) + a0) - b0)
# If clamping: the only time we'll get closest points will be when lines don't overlap at all.
# Find if segments overlap using dot products.
if clampA0 or clampA1 or clampB0 or clampB1:
d1 = np.dot(_A, (b1 - a0))
# Is segment B before A?
if d0 <= 0 >= d1:
if clampA0 == True and clampB1 == True:
if np.absolute(d0) < np.absolute(d1):
return b0, a0, np.linalg.norm(b0 - a0)
return b1, a0, np.linalg.norm(b1 - a0)
# Is segment B after A?
elif d0 >= np.linalg.norm(A) <= d1:
if clampA1 == True and clampB0 == True:
if np.absolute(d0) < np.absolute(d1):
return b0, a1, np.linalg.norm(b0 - a1)
return b1, a1, np.linalg.norm(b1, a1)
# If clamping is off, or segments overlapped, we have infinite results, just return position.
return None, None, d
# Lines criss-cross: Calculate the dereminent and return points
t = (b0 - a0);
det0 = np.linalg.det([t, _B, cross])
det1 = np.linalg.det([t, _A, cross])
t0 = det0 / denom;
t1 = det1 / denom;
pA = a0 + (_A * t0);
pB = b0 + (_B * t1);
# Clamp results to line segments if needed
if clampA0 or clampA1 or clampB0 or clampB1:
if t0 < 0 and clampA0:
pA = a0
elif t0 > np.linalg.norm(A) and clampA1:
pA = a1
if t1 < 0 and clampB0:
pB = b0
elif t1 > np.linalg.norm(B) and clampB1:
pB = b1
d = np.linalg.norm(pA - pB)
return pA, pB, d
'''
Calculate the distance of the line segment between the start and end points
'''
def distanceBetweenTwoPoints(startPoint, endPoint):
return np.linalg.norm(endPoint - startPoint)
'''
A function used to reduce a line segment between the startPoint and the endPoint
by the recduce amount. So if you have a line that is 1 distance long and you want to
move both of its end points in by .25 this function will return the two new end points
and the length of the new line will be .5.
'''
def reduceLineSegment(startPoint, endPoint, reduceAmount):
vector = []
'''
Parametrize the line segment so that you can subtract the distace. Since the line
is parametrized its distance will be from 0 to 1. To reduce it by the desired amount
we need to figure out the proportion of the line compared to the parameter and then
subtract that amount.
'''
vector.append(startPoint[0] - endPoint[0]) # Calculate Vector X
vector.append(startPoint[1] - endPoint[1]) # Calculate Vector Y
vector.append(startPoint[2] - endPoint[2]) # Calculate Vector Z
dist = distanceBetweenTwoPoints(startPoint, endPoint)
proportion = reduceAmount / dist
newX = endPoint[0] + proportion * vector[0]
newY = endPoint[1] + proportion * vector[1]
newZ = endPoint[2] + proportion * vector[2]
newEndPoint = np.array([newX, newY, newZ])
output = []
newX = endPoint[0] + (1 - proportion) * vector[0]
newY = endPoint[1] + (1 - proportion) * vector[1]
newZ = endPoint[2] + (1 - proportion) * vector[2]
newStartPoint = np.array([newX, newY, newZ])
output.append(newStartPoint)
output.append(newEndPoint)
return output
'''
Generates the set of second iterated forward difference operator. This comes from Professor Peters' and
Ji Li's paper. In the paper it is define. It shows up as a triangle with a subscript 2 and and set P next to it
'''
def secondIteratedForwardDifferenceOperator(listOfControlPoints, index):
deltaSub2Set = []
numCtrlPts = len(listOfControlPoints)
# Need to wrap around and connect back to the first point
for i in range(0, numCtrlPts): # Add -1 to numCtrlPts to make it an open curve again
deltaSub2Ofi = listOfControlPoints[(i + 2) % numCtrlPts][index] - (2 * listOfControlPoints[(i + 1) % numCtrlPts][index]) + listOfControlPoints[i][index]
deltaSub2Set.append(deltaSub2Ofi)
return deltaSub2Set
'''
This is the l1 - Norm for a set. The l1 norm is just the summation of the absolute value of all the
elements in the set.
'''
def l1Norm(setForNorming):
norm = 0
i = 0
for element in setForNorming:
norm += abs(element)
i += 1
print('Terms in L1 Norm: ', i)
return norm
'''
This generates the Omega one value from Professor Peters' and Ji Li's paper. Omega one is used in combination with the delta from the Denne-Sullivan paper to figure out m1
'''
def omegaOne(listOfControlPoints):
deltaSetX = secondIteratedForwardDifferenceOperator(listOfControlPoints, 0)
deltaSetY = secondIteratedForwardDifferenceOperator(listOfControlPoints, 1)
deltaSetZ = secondIteratedForwardDifferenceOperator(listOfControlPoints, 2)
l1NormX = l1Norm(deltaSetX)
l1NormY = l1Norm(deltaSetY)
l1NormZ = l1Norm(deltaSetZ)
return max(l1NormX, l1NormY, l1NormZ)
'''
This function generates the m1 value from Professor Peters' and Ji Li's paper. M1 is compared against
m2 and m3 to determine the number of iterations, j, that need to be done in order for the stick knot
to properly associate with its correct bezier curve
'''
def generateM1(omegaOne, delta):
# omegaOneSquared = omegaOne * omegaOne
omegaOneSquared = np.multiply(omegaOne, omegaOne)
deltaSquared = delta * delta
intermediate = ((float(7) / float(16) * omegaOneSquared * (1 / deltaSquared)) - (float(1) / float(7)))
logResult = math.log(intermediate, 2)
# return math.ceil(logResult)
intermediate = np.log2((7 * omegaOne) - (4 * delta)) + np.log2((7 * omegaOne) + (4 * delta)) - np.log2(112 * delta * delta)
answer = np.ceil(intermediate)
# TODO: PDZ- This needs to NOT be hardcoded. Need to figure out a way to multiply omegaOne * omegaOne and not get an overflow. Then the two lines above can be uncommented.
# return math.ceil(90.958811263220) # Unknot hard coded
return math.ceil(104.34285901943) # Trefoil Hard coded
'''
Generates the hodograph delta sub 2 set from Professor Peters' and Ji Li's paper. It is defined in the paper
'''
def generateHodographDeltaTwoSet(listOfControlPoints, index):
hodographDelta2Set = []
numControlPoints = len(listOfControlPoints)
# Need to wrap around and connect back to the first point
for i in range(1, numControlPoints): # Add -1 to numControlPoints to make it an open curve again
hodographDeltaSub2Ofi = numControlPoints * (listOfControlPoints[(i + 2) % numControlPoints][index] - (
3 * listOfControlPoints[(i + 1) % numControlPoints][index]) + (3 * listOfControlPoints[i][index]) -
listOfControlPoints[i - 1][index])
hodographDelta2Set.append(hodographDeltaSub2Ofi)
return hodographDelta2Set
'''
Found the same way as Omega One but uses the hodograph sets instead
'''
def omegaTwo(listOfControlPoints):
hodographX = generateHodographDeltaTwoSet(listOfControlPoints, 0)
hodographY = generateHodographDeltaTwoSet(listOfControlPoints, 1)
hodographZ = generateHodographDeltaTwoSet(listOfControlPoints, 2)
l1NormX = l1Norm(hodographX)
l1NormY = l1Norm(hodographY)
l1NormZ = l1Norm(hodographZ)
return max(l1NormX, l1NormY, l1NormZ)
'''
The M2 value from Professor Peters' and Ji Li's paper. It is found using the hodograph set
'''
def generateM2(omegaTwo, lambdaVal, numCtrlPts):
j = 0
# This gets the next int higher to what would be an unreducible logarithm
while ((numCtrlPts * math.pow(2, 3 * j) + math.pow(2, 2 * j)) < math.pow(omegaTwo / lambdaVal, 2)):
j += 1
return j
'''
The M3 value from Professor Peters' and Ji Li's paper. It is very similar to M2 except for the additional sin(pi/ 8) value
'''
def generateM3(omegaTwo, lambdaVal, numCtrlPts):
j = 0
# This gets the next int higher to what would be an unreducible logarithm
while ((numCtrlPts * math.pow(2, 3 * j) + math.pow(2, 2 * j)) < math.pow(
(omegaTwo / (math.sin(math.pi / 8) * lambdaVal)), 2)):
j += 1
return j
if __name__ == '__main__':
np.seterr(all='raise', over='raise')
[testa, testb, testdist] = closestDistanceBetweenLines(np.array([0.0, 0.0, 0.0]), np.array([5.0, 0.0, 0.0]),
np.array([6.0, 0.0, 1.0]), np.array([10.0, 0.0, 1.0]),
clampAll=True)
print("Done: ", testdist)
'''
Add the line segments of the control polygon to a list
'''
# The Unknot
segmentArray = []
# segmentArray.append(np.array([0.0, 4831838208, 10737418240]))
# segmentArray.append(np.array([-8053063680, -51002736640, -26843545600]))
# segmentArray.append(np.array([21474836480, 42949672960, -10737418240]))
# segmentArray.append(np.array([-5368709120, -32212254720, 31138512896]))
# segmentArray.append(np.array([-32212254720, 16106127360, 10737418240]))
# segmentArray.append(np.array([21474836480, -32212254720, -32212254720]))
# #segmentArray.append(np.array([0, 4831838208,10737418240]))
# The Trefoil
segmentArray.append(np.array([0.0, 4831838208, 10737418240]))
segmentArray.append(np.array([-8053063680, -51002736640, -26843545600]))
segmentArray.append(np.array([21474836480, 42949672960, -10737418240]))
segmentArray.append(np.array([5368709120, -32212254720, 31138512896]))
segmentArray.append(np.array([-32212254720, 16106127360, 10737418240]))
segmentArray.append(np.array([21474836480, -32212254720, -32212254720]))
[a, b, dist] = closestDistanceBetweenLines(segmentArray[0], segmentArray[1], segmentArray[2], segmentArray[3],
clampAll=True)
print('TEST: a is: ', a, ' b is: ', b, ' dist is: ', dist)
'''
segmentArray = []
segmentArray.append(np.array([1.28318923, -2.84316645, -2.25593748]))
segmentArray.append(np.array([1.139367593750000, -2.831090593750000, -2.293403687500000]))
segmentArray.append(np.array([0.971135187500000, -2.330181187500000, -2.079607375000000]))
segmentArray.append(np.array([0.802902781250000, -1.829271781250000, -1.865811062500000]))
segmentArray.append(np.array([0.634670375000000, -1.328362375000000, -1.652014750000000]))
segmentArray.append(np.array([0.466437968750000, -0.827452968750000, -1.438218437500000]))
segmentArray.append(np.array([0.298205562500000, -0.326543562500000, -1.224422125000000]))
segmentArray.append(np.array([0.129973156250000, 0.174365843750000, -1.010625812500000]))
segmentArray.append(np.array([-0.038259250000000, 0.675275250000000, -0.796829500000000]))
segmentArray.append(np.array([-0.206491656250000, 1.176184656250000, -0.583033187500000]))
segmentArray.append(np.array([-0.374724062500000, 1.677094062500000, -0.369236875000000]))
segmentArray.append(np.array([-0.542956468750000, 2.178003468750000, -0.155440562500000]))
segmentArray.append(np.array([-0.711188875000000, 2.678912875000000, 0.058355750000000]))
segmentArray.append(np.array([-0.879421281250000, 3.179822281250000, 0.272152062500000]))
segmentArray.append(np.array([-1.047653687500000, 3.680731687500000, 0.485948375000000]))
segmentArray.append(np.array([-1.215886093750000, 4.181641093750000, 0.699744687500000]))
segmentArray.append(np.array([-1.384118500000000, 4.682550500000000, 0.913541000000000]))
segmentArray.append(np.array([-1.50375, 4.13635, 1.02434]))
segmentArray.append(np.array([-1.62339, 3.59015, 1.13513]))
segmentArray.append(np.array([-1.74303, 3.04394, 1.24592]))
segmentArray.append(np.array([-1.86266, 2.49774, 1.35671]))
segmentArray.append(np.array([-1.9823, 1.95154, 1.4675]))
segmentArray.append(np.array([-2.10194, 1.40534, 1.57829]))
segmentArray.append(np.array([-2.22157, 0.859137, 1.68908]))
segmentArray.append(np.array([-2.34121, 0.312934, 1.79987]))
segmentArray.append(np.array([-2.46084, -0.233267, 1.91066]))
segmentArray.append(np.array([-2.58048, -0.779468, 2.02145]))
segmentArray.append(np.array([-2.70012, -1.32567, 2.13224]))
segmentArray.append(np.array([-2.81976, -1.87187, 2.24303]))
segmentArray.append(np.array([-2.9394, -2.41807, 2.35382]))
segmentArray.append(np.array([-3.05903, -2.96428, 2.46461]))
segmentArray.append(np.array([-3.17867, -3.51048, 2.5754]))
segmentArray.append(np.array([-3.2983075,-4.0566825,2.686189]))
segmentArray.append(np.array([-3.09987, -3.63013, 2.36433]))
segmentArray.append(np.array([-2.90143, -3.20357, 2.04247]))
segmentArray.append(np.array([-2.70299, -2.77701, 1.72061]))
segmentArray.append(np.array([-2.50455, -2.35045, 1.39875]))
segmentArray.append(np.array([-2.30611, -1.92389, 1.07689]))
segmentArray.append(np.array([-2.10767, -1.49733, 0.755026]))
segmentArray.append(np.array([-1.90924, -1.07077, 0.433165]))
segmentArray.append(np.array([-1.7108, -0.644214, 0.111303]))
segmentArray.append(np.array([-1.51236, -0.217655, -0.210558]))
segmentArray.append(np.array([-1.31393, 0.208903, -0.532419]))
segmentArray.append(np.array([-1.11549, 0.635462, -0.854279]))
segmentArray.append(np.array([-0.91705, 1.06202, -1.17614]))
segmentArray.append(np.array([-0.718613, 1.48858, -1.498]))
segmentArray.append(np.array([-0.520175, 1.91514, -1.81986]))
segmentArray.append(np.array([-0.321737, 2.3417, -2.14172]))
segmentArray.append(np.array([-0.1232995,2.768254,-2.463584]))
segmentArray.append(np.array([0.128657, 2.3119, -2.23296]))
segmentArray.append(np.array([0.380613, 1.85555, -2.00234]))
segmentArray.append(np.array([0.632569, 1.3992, -1.77172]))
segmentArray.append(np.array([0.884525, 0.942852, -1.5411]))
segmentArray.append(np.array([1.13648, 0.486501, -1.31047]))
segmentArray.append(np.array([1.38844, 0.0301505, -1.07985]))
segmentArray.append(np.array([1.64039, -0.4262, -0.849228]))
segmentArray.append(np.array([1.89235, -0.882551, -0.618607]))
segmentArray.append(np.array([2.14431, -1.3389, -0.387985]))
segmentArray.append(np.array([2.39626, -1.79525, -0.157363]))
segmentArray.append(np.array([2.64821, -2.2516, 0.0732595]))
segmentArray.append(np.array([2.90017, -2.70795, 0.303882]))
segmentArray.append(np.array([3.15212, -3.1643, 0.534504]))
segmentArray.append(np.array([3.40408, -3.62065, 0.765126]))
segmentArray.append(np.array([3.65604, -4.077, 0.995748]))
segmentArray.append(np.array([3.9079915,-4.533357,1.2263705]))
segmentArray.append(np.array([3.41775, -4.27741, 1.08826]))
segmentArray.append(np.array([2.9275, -4.02147, 0.950154]))
segmentArray.append(np.array([2.43725, -3.76553, 0.812045]))
segmentArray.append(np.array([1.947, -3.50958, 0.673937]))
segmentArray.append(np.array([1.45675, -3.25364, 0.535829]))
segmentArray.append(np.array([0.966502, -2.9977, 0.39772]))
segmentArray.append(np.array([0.476253, -2.74175, 0.259611]))
segmentArray.append(np.array([-0.0139957, -2.48581, 0.121503]))
segmentArray.append(np.array([-0.504244, -2.22986, -0.0166055]))
segmentArray.append(np.array([-0.994493, -1.97392, -0.154714]))
segmentArray.append(np.array([-1.48474, -1.71798, -0.292822]))
segmentArray.append(np.array([-1.97499, -1.46204, -0.430931]))
segmentArray.append(np.array([-2.46524, -1.2061, -0.56904]))
segmentArray.append(np.array([-2.95549, -0.950156, -0.707148]))
segmentArray.append(np.array([-3.44574, -0.694214, -0.845257]))
segmentArray.append(np.array([-3.935983,-0.438272,-0.983365]))
segmentArray.append(np.array([-3.48885, -0.142371, -0.789876]))
segmentArray.append(np.array([-3.04171, 0.153529, -0.596387]))
segmentArray.append(np.array([-2.59457, 0.449429, -0.402898]))
segmentArray.append(np.array([-2.14744, 0.745329, -0.209409]))
segmentArray.append(np.array([-1.7003, 1.04123, -0.01592]))
segmentArray.append(np.array([-1.25317, 1.33713, 0.177569]))
segmentArray.append(np.array([-0.806037, 1.63303, 0.371058]))
segmentArray.append(np.array([-0.358904, 1.92893, 0.564547]))
segmentArray.append(np.array([0.0882295, 2.22483, 0.758035]))
segmentArray.append(np.array([0.535363, 2.52073, 0.951523]))
segmentArray.append(np.array([0.982496, 2.81663, 1.14501]))
segmentArray.append(np.array([1.42963, 3.11253, 1.3385]))
segmentArray.append(np.array([1.87677, 3.40843, 1.53199]))
segmentArray.append(np.array([2.3239, 3.70433, 1.72548]))
segmentArray.append(np.array([2.77104, 4.00023, 1.91897]))
segmentArray.append(np.array([3.218174000000000, 4.296123000000000, 2.112459500000000]))
segmentArray.append(np.array([3.098763125000000, 3.819365312500000, 1.823730781250000]))
segmentArray.append(np.array([2.979352250000000, 3.342607625000000, 1.535002062500000]))
segmentArray.append(np.array([2.859941375000000, 2.865849937500000, 1.246273343750000]))
segmentArray.append(np.array([2.740530500000000, 2.389092250000000, 0.957544625000000]))
segmentArray.append(np.array([2.621119625000000, 1.912334562500000, 0.668815906250000]))
segmentArray.append(np.array([2.501708750000000, 1.435576875000000, 0.380087187500000]))
segmentArray.append(np.array([2.382297875000000, 0.958819187500000, 0.091358468750000]))
segmentArray.append(np.array([2.262887000000000, 0.482061500000000, -0.197370250000000]))
segmentArray.append(np.array([2.143476125000000, 0.005303812500000, -0.486098968750000]))
segmentArray.append(np.array([2.024065250000000, -0.471453875000000, -0.774827687500000]))
segmentArray.append(np.array([1.904654375000000, -0.948211562500000, -1.063556406250000]))
segmentArray.append(np.array([1.785243500000000, -1.424969250000000, -1.352285125000000]))
segmentArray.append(np.array([1.665832625000000, -1.901726937500000, -1.641013843750000]))
segmentArray.append(np.array([1.546421750000000, -2.378484625000000, -1.929742562500000]))
segmentArray.append(np.array([1.427010875000000, -2.855242312500000, -2.218471281250000]))
'''
minimums = []
'''
For each line segment, compare it to all non-incident (non-adjacent) line segments. Add the distances
to the minimums list.
'''
numCtrlPointSegments = len(segmentArray)
for i in range(0, numCtrlPointSegments):
# print'i is: ' ,i
j = i + 2
while j < numCtrlPointSegments:
if i == 0 and (j + 1) % numCtrlPointSegments == 0:
j += 1
continue
# print 'j is: ', j, ' and j + 1 mod 7 is', (j+1)% numCtrlPointSegments
[a, b, dist] = closestDistanceBetweenLines(segmentArray[i], segmentArray[i + 1], segmentArray[j],
segmentArray[(j + 1) % numCtrlPointSegments], clampAll=True)
print('a is: ', a, ' b is: ', b, ' dist is: ', dist)
# print(dist)
minimums.append(dist)
j += 1
'''
Of all the distances, take the minimum. This is the r1 value from the Denne-Sullivan paper
'''
r1 = min(minimums)
print('r1 is: ', r1)
'''
Epsilon is a constant that we have pre selected and can change
'''
epsilon = 1.0
'''
r2 is also from the Denne-Sullivan paper
'''
r2 = min(r1 / 2, epsilon / 2)
print('r2 is: ', r2)
'''
For each line segment of the control polygon, remove r2 from each end.
'''
adjustedSegments = []
for i in range(0, numCtrlPointSegments):
# print 'i is: ', i, 'and the Line Segment is ', segmentArray[i], segmentArray[(i+1)%numCtrlPointSegments]
[segment1, segment2] = reduceLineSegment(segmentArray[i], segmentArray[(i + 1) % numCtrlPointSegments], r2)
adjustedSegments.append(segment1)
adjustedSegments.append(segment2)
# print 'The new segment is: ', segment1, segment2
i += 1
# At this point the segments have been reduced along the vertex. This means that the number of points is 2 times the number of control points. This is because the segments no longer share any common points between themselves
numAdjustedPoints = len(adjustedSegments)
newMinimums = []
'''
Then, for each reduced segment find the distance from that segment to all other reduced segments
'''
i = 0
while i + 3 < numAdjustedPoints:
j = i + 2
while j + 1 < numAdjustedPoints:
# print 'Comparing: ', adjustedSegments[i], adjustedSegments[i+1], 'to: ', adjustedSegments[j],adjustedSegments[j+1]
[a, b, dist] = closestDistanceBetweenLines(adjustedSegments[i], adjustedSegments[i + 1],
adjustedSegments[j], adjustedSegments[j + 1], clampAll=True)
newMinimums.append(dist)
# print dist
j += 2
i += 2
'''
Take r3 (also from Denne-Sullivan) to be the minimum distance of all distances between line segments
'''
r3 = min(newMinimums)
print('r3 is: ', r3)
'''
r4 is also from Denne-Sullivan paper
'''
r4 = r3 / 6
print('r4 is: ', r4)
'''
Finally, calculate delta (from Denne-Sullivan)
'''
delta = r4 / 3
print('delta is: ', delta)
omega1 = omegaOne(segmentArray)
print('Omega One is: ', omega1)
m1 = generateM1(omega1, delta)
print('M1 is: ', m1)
segmentLengths = []
for i in range(0, numCtrlPointSegments):
segmentLengths.append(distanceBetweenTwoPoints(segmentArray[i], segmentArray[(i + 1) % numCtrlPointSegments]))
lambdaVal = min(segmentLengths)
print('Lambda is: ', lambdaVal)
omega2 = omegaTwo(segmentArray)
print('Omega Two is: ', omega2)
m2 = generateM2(omega2, lambdaVal, numCtrlPointSegments)
print('M2 is: ', m2)
m3 = generateM3(omega2, lambdaVal, numCtrlPointSegments)
print('M3 is: ', m3)