0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022"""inherit from cvisual vector class, and add geometric methods, serving
0023the analytics for real points in space"""
0024
0025
0026from pygeo.base.abstracts._element import Element
0027from pygeo.base.analytics.pygeomath import *
0028from pygeo.base.support.pygeoopts import EPS, MAX ,DO_TESTS
0029
0030class Position3(vector):
0031 """a 3 element vector representing a position in 3d space"""
0032
0033 def __init__(self,*args):
0034 """filter out the non-numeric arguments of a derived class, and
0035 send the numeric arguments to the initialization of the vector parent
0036 class"""
0037
0038 v=[]
0039 for arg in args:
0040 if not isinstance(arg,Element):
0041 v.append(arg)
0042 else:
0043 if isinstance(arg,vector):
0044 v=vector(arg.x,arg.y,arg.z)
0045
0046
0047
0048 vector.__init__(self,*v)
0049
0050 def set(self,other):
0051 """ set defining attributes of point to a given point"""
0052
0053 self.x=other.x
0054 self.y=other.y
0055 self.z=other.z
0056 def polar(self):
0057 """ return the polar coordinates of the point"""
0058
0059 s = hypot(self.x,self.y)
0060 t=arctan2(self.y,self.x)
0061 p=arctan2(s,self.z)
0062 return p,t
0063
0064 def distance(self,other):
0065 """ the distance between the point and a second point"""
0066
0067 return (self-other).mag
0068
0069 def distanceSquared(self,other):
0070 """ the square of the distance between the point and a second point"""
0071
0072 return (self-other).mag2
0073
0074 def homogenous(self):
0075 """ a flat array representing the homogenous coordinates of the point"""
0076
0077 return array((self.x,self.y,self.z,1.),'d')
0078
0079 def transform(self,mat,point):
0080 """return the given point as 'self' transformed by a given a 4x4 transformation matrix"""
0081
0082 tf=matrixmultiply(homogenous(self),mat)
0083 point.set(vector(tf[0:3]/tf[3]))
0084
0085 def onPlane(self,plane):
0086 """ is point on a given plane?"""
0087
0088 return absolute(matrixmultiply(plane.equat(),self.homogenous())) < EPS
0089
0090 def coPlanar(self,p1,p2,p3):
0091 """is point coplanar with 3 given points?"""
0092
0093 c= absolute(determinant(array([p1.homogenous(),
0094 p2.homogenous(),
0095 p3.homogenous(),
0096 self.homogenous()])))
0097
0098 if c < EPS:
0099 return True
0100 else:
0101 return False
0102
0103 def coLinear(self,p1,p2):
0104 """ is point colinear with 2 given points ?"""
0105
0106 v1=p1-self
0107 v2=p2-p1
0108 return v1.cross(v2).mag2 < EPS
0109
0110 def onLine(self,line):
0111 """ is point on given line?"""
0112
0113 v1=line.p1 - self
0114 return v1.cross(line.getDirection).mag2 < EPS
0115
0116 def toInterpolated(self,p1,p2,r):
0117 """ set the position on a line connecting to given points,
0118 based on the interpolation by a given ratio of the distance between
0119 the given points"""
0120
0121 self.set(p1*(1.0-r)+p2*r)
0122 return True
0123
0124 def to_3d(self,harray):
0125 """ set position as the projection from given homogenous coordinates"""
0126
0127 try :
0128 self.set(vector(harray[0:3]/float(harray[3])))
0129 except (ZeroDivisionError,OverflowError):
0130 print self.__class__.__name__
0131 print "point at infinity in conversion from homogenous coordinates"
0132 print "vector set to 'MAX' times non-homogenous coordinates"
0133 self.set(vector(harray[0:3]*MAX))
0134 return True
0135
0136 def toXY(self,point):
0137 """ set position as the origin centered projection of the given point
0138 to the plane parallel with the XY plane at unit distance"""
0139
0140 __mat2=array([[-1.,0.,0.,0.],[0.,-1.,0.,0.],[0.,0.,0.,-1.],[0.,0.,0.,-1]])
0141 tf= matrixmultiply(point.homogenous(),__mat2)
0142 self.set(vector([tf[0]/tf[3],tf[1]/tf[3],1]))
0143 return True
0144
0145 def fromXY(self,equat,point):
0146 """ set position as the origin centered projection of the given point
0147 to the plane represented by the given 4 element array"""
0148
0149 _zpt=array([0.,0.,-1.,1.])
0150 mat= multiply.outer(equat,_zpt)
0151 k=matrixmultiply(_zpt,equat)
0152 for i in range(4):
0153 mat[i,i]-=k
0154 tf=matrixmultiply(homogenous(point),mat)
0155 try :
0156 self.set(vector(tf[0:3]/tf[3]))
0157 except (ZeroDivisionError,OverflowError):
0158 print self.__class__.__name__
0159 print "point at infinity in conversion from homogenous coordinates"
0160 print "vector set to 'MAX' times non-homogenous coordinates"
0161 self.set(vector(tf[0:3][0:3])*MAX)
0162 return True
0163
0164 def toLine(self,p1,p2):
0165 """move to foot of perpendicular from initial position of self
0166 to line connecting p1 & p2 arguments"""
0167
0168 vxs = p2 - p1
0169 vxt = self - p1
0170 try:
0171 factor= (
0172 vxt.dot(vxs)
0173 /
0174 vxs.mag2
0175 )
0176 except ZeroDivisionError:
0177 print self.__class__.__name__
0178 print " points defining line segments are coincident. toLine returned False"
0179 return False
0180 self.set(vxs*factor+p1)
0181 return True
0182
0183 def toPlane(self,plane):
0184 """move to foot of perpendicular from initial position of self
0185 to plane argument"""
0186
0187 pt=homogenous(self-plane._u)
0188 equat=plane.equat()
0189 mat= multiply.outer(equat,pt)
0190 k=matrixmultiply(pt,equat)
0191 for i in range(4):
0192 mat[i,i]-=k
0193 self.to_3d(matrixmultiply(self.homogenous(),mat))
0194 return True
0195
0196 def toCircle(self,circle):
0197 """move to point on the given circle's circumference closest to the foot of
0198 perpendicular from initial position of self to the circle's plane"""
0199
0200 self.toPlane(circle)
0201 vt=self-circle._center
0202 if vt.mag:
0203 factor = circle._radius/vt.mag
0204 self.set(vt*factor+circle._center)
0205 else:
0206 lxy = hypot(circle._u.x,circle._u.y)
0207 if (lxy >= EPS):
0208 self.set(vector(-circle._u.y/lxy,circle._u.x/lxy,0).norm()
0209 *circle._radius)
0210 else:
0211 self.vector=vector(1,0,0)
0212 return True
0213
0214 def toSphere (self,sphere):
0215 """move to closest point on surface of sphere argument"""
0216
0217 vt= self- sphere._center
0218 len=vt.mag
0219 if len:
0220 factor=sphere._radius/len
0221 self.set(vt*factor+sphere._center)
0222 else:
0223 self.set(vector(sphere._center.x+sphere._radius,
0224 sphere._center.y,
0225 sphere._center.z))
0226 return True
0227
0228 def toCentroid(self,p1,p2,p3):
0229 """ move to the centroid of the triangle defined by the 3 point arguments"""
0230
0231 self.set((p1+p2+p3)/3.0)
0232 return True
0233
0234 def toOrthoCenter(self,p1,p2,p3):
0235 """ move to the orthocenter of the triangle defined by the 3 point arguments"""
0236
0237 vxs=Position3()
0238 vxs.set(p3)
0239 vxs.toLine(p2,p1)
0240 vxt=Position3()
0241 vxt.set(p2)
0242 vxt.toLine(p3,p1)
0243 if self.toInterSection(p3,vxs,p2,vxt,test=False):
0244 return True
0245 else:
0246 self.__class__.__name__
0247 print "points are not distinct, no orthocenter defined, returned False"
0248 return False
0249
0250 def toInCenter(self,p1,p2,p3):
0251 """ move to the center of the circle inscribed in the
0252 triangle defined by the 3 point arguments"""
0253
0254 vxs=Position3()
0255 vxt=Position3()
0256 vxs.toBisector(p1,p2,p3)
0257 vxt.toBisector(p2,p1,p3)
0258 if self.toInterSection(p1,vxs,p2,vxt,test=False):
0259 return True
0260 else:
0261 self.__class__.__name__
0262 print "points are not distinct, no incenter defined, returned False"
0263 return False
0264
0265 def toExCenter(self,p1,p2,p3):
0266 """ move to the excenter of the circle exscribed in the
0267 triangle defined by the 3 point arguments, and tangent to the side
0268 opposite the first point argument"""
0269
0270 vxs=Position3()
0271 vxt=Position3()
0272 vxs.toExtBisector(p2,p1,p3)
0273 vxt.toExtBisector(p3,p1,p2)
0274 if self.toInterSection(p2,vxs,p3,vxt,test=0):
0275 return True
0276 else:
0277 self.__class__.__name__
0278 print "points are not distinct, no excenter defined, returned False"
0279 return False
0280
0281 def toOpposite(self,object,point):
0282 """ set postion as antipodal to the given point relative to the given object, i.e.
0283 circle or sphere"""
0284
0285 try:
0286 v= (object._center-point).norm()
0287 self.set(v*object._radius+object._center)
0288 return True
0289 except ZeroDivisionError:
0290 print self.__class__.__name__
0291 print "point at object center, antipodal undefined, returned False"
0292 return False
0293
0294 def toInterSection(self,p11,p12,p21,p22,test=True):
0295 """ set position as the intersection of the line connecting the
0296 1st and 2nd point arguments with that connecting the 3rd and 4th
0297 point arugments. Default is to test if lines are slew"""
0298
0299 if DO_TESTS:
0300 if test:
0301 if not p11.coPlanar(p12,p21,p22):
0302 print self.__class__.__name__
0303 print "no intersection of skew lines; False returned"
0304 return False
0305 d1 = p12 - p11
0306 d2 = p22 - p21
0307 c = d1.cross(d2)
0308 o = p21 - p11
0309 try:
0310 factor = determinant(array((o,d2,c)))/c.mag2
0311 except ZeroDivisionError:
0312 if (d1.mag <> 0 and d2.mag <> 0):
0313 print self.__class__.__name__
0314 print "intersection of parallel lines"
0315 print "returning MAX times vector"
0316 factor = MAX
0317 else:
0318 return False
0319 self.set(d1*factor+p11)
0320 return True
0321
0322 def toPlaneIntersection(self,plane,p1,p2,test=False):
0323 """set to position of intersection of given plane with the
0324 line connecting the point arguments. Default is not to test if line
0325 is on plane"""
0326
0327 if test:
0328 if p1.onPlane(plane) and p2.onPlane(plane):
0329 print self.__class__.__name__
0330 print "line on plane, degenerate intersection, False returned"
0331 return False
0332 pt1=p1.homogenous()
0333 pt2=p2.homogenous()
0334 equat=plane.equat()
0335 mat= multiply.outer(equat,pt2)
0336 k=matrixmultiply(pt2,equat)
0337 for i in range(4):
0338 mat[i,i]-=k
0339 self.to_3d(matrixmultiply(pt1,mat))
0340 return True
0341
0342 def toPlaneReflection(self,u,d):
0343 """move to the reflection of the current position in
0344 the plane defined by the 3 element vector representing the
0345 planes unit normal and a numeric value representing the plane's
0346 distance from the origin"""
0347
0348 x,y,z=u.x,u.y,u.z
0349 mat=transpose(array([[1-2*x**2,-2*x*y,-2*x*z,2*x*d],
0350 [-2*x*y,1-2*y**2,-2*y*z,2*y*d],
0351 [-2*x*z,-2*y*z,1-2*z**2,2*z*d],
0352 [0.,0.,0.,1]]))
0353 self.to_3d(matrixmultiply(self.homogenous(),mat))
0354 return True
0355
0356 def toSphereCenter(self,p1,p2,p3,p4):
0357 """ set position as the center of the sphere defined by the 4 given
0358 points"""
0359
0360 s1=p1.mag2
0361 s2=p2.mag2
0362 s3=p3.mag2
0363 s4=p4.mag2
0364 v21=p2-p1
0365 v31=p3-p1
0366 v41=p4-p1
0367 n1=determinant(array(((s2-s1,v21.y,v21.z),(s3-s1,v31.y,v31.z),(s4-s1,v41.y,v41.z))))
0368 n2=determinant(array(((v21.x,s2-s1,v21.z),(v31.x,s3-s1,v31.z),(v41.x,s4-s1,v41.z))))
0369 n3=determinant(array(((v21.x,v21.y,s2-s1),(v31.x,v31.y,s3-s1),(v41.x,v41.y,s4-s1))))
0370 d3=determinant(array((v21,v31,v41)))
0371 try:
0372 self.set(vector(float(n1)/(d3*2),float(n2)/(d3*2),float(n3)/(d3*2)))
0373 return True
0374 except ZeroDivisionError:
0375 print self.__class__.__name__
0376 print "points are not distinct, sphere center undefined, returned False"
0377 return False
0378
0379 def toCircumCenter(self,p1,p2,p3):
0380 """ set position as the center of the circle on the 3 given points"""
0381
0382 vxs = p2 - p1
0383 vxt = p3 - p1
0384 d1=vxs.dot(vxt)
0385 d2=vxs.mag2
0386 d3= vxt.mag2
0387 den=2.0*(d2*d3-d1*d1)
0388 try:
0389 f1=(d2-d1)/den*d3
0390 f2=(d3-d1)/den*d2
0391 except ZeroDivisionError:
0392 print self.__class__.__name__
0393 print "points are not distinct, circumcenter not defined, returned False"
0394 return False
0395 self.set(vxs*f1+vxt*f2+p1)
0396 return True
0397
0398 def toBisector(self,p1,p2,p3):
0399 """set position the (interior) bisector of the angle at the 1st point argument
0400 formed with the the 2nd and 3rd point arguments"""
0401 try:
0402 r=1.0/(p1.distance(p3)/p1.distance(p2)+ 1.0)
0403 except ZeroDivisionError:
0404 r=1.0
0405 self.toInterpolated(p2,p3,r)
0406 return True
0407
0408 def toExtBisector(self,p1,p2,p3):
0409 """set position the exterior bisector of the angle at the 1st point argument
0410 formed with the the 2nd and 3rd point arguments"""
0411
0412 try:
0413 r=1.0/(-p1.distance(p3)/p1.distance(p2)+1.0)
0414 except ZeroDivisionError:
0415 r=1.0
0416 self.toInterpolated(p2,p3,r)
0417 return True
0418
0419 def toInvertPoint(self,circle,p1):
0420 """set position as the inversion of the given point in the
0421 given circle, on the circle's plane"""
0422
0423 t=(p1 - circle._center)*circle._radiusSquared
0424 try:
0425 t /= float(p1.distanceSquared(circle._center))
0426 self.set(t + circle._center)
0427 except ZeroDivisionError:
0428 print self.__class__.__name__
0429 print "point at infinity returned by 'toInvertPoint', returning MAX"
0430 self.set(vector(MAX,MAX,MAX))
0431 return True
0432
0433 def toCrossPoint(self,p1,p2,p3,p4,p1a,p2a,p3a):
0434 """set position as the crosspoint (equating cross ratios) of the
0435 last 3 (ordered) points with the first 4 (ordered) points"""
0436
0437 v1=vector(multiply((p3-p1),(p4-p2)))
0438 d1=v1.dot(p3a-p2a)
0439 v2=vector(multiply((p3-p2),(p4-p1)))
0440 d2=v2.dot(p3a-p1a)
0441 try:
0442 d=d1/d2
0443 except ZeroDivisionError:
0444 d=MAX
0445 try:
0446 t = 1.0/(1.0 - d)
0447 except ZeroDivisionError:
0448 t=0
0449 self.toInterpolated(p1a,p2a,t)
0450 return True
0451
0452 def toCircumPoint(self,circle,rad):
0453 """ rotate the current position (assumed to be on the circle's
0454 circumference) by the given angle, in radians"""
0455
0456 if not self==circle._center:
0457 self.set((self-circle._center).rotate(rad,
0458 circle._u) + circle._center)
0459 return True
0460 else:
0461 return False
0462
0463 def toHarmonic(self,p1,p2,p3):
0464 """ set position as the harmonic with of the 3rd point argument with reespect
0465 to the 1st and 2nd point arguments"""
0466
0467 t3=(p2 + p3) * .5
0468 t2= p2 - t3
0469 t1 = p1 -t3
0470 m1 = t1.mag2
0471 m2 = t2.mag2
0472 v=t2*(2.0*t1.dot(t2)) - t1*m2
0473 try:
0474 self.set(v/m1+t3)
0475 except ZeroDivisionError:
0476 print self.__class__.__name__
0477 print "harmonic to mipoint, point at infinity, returned [MAX,MAX,MAX}"
0478 self.set(vector(MAX,MAX,MAX))
0479 return True
0480
0481
0482__author__ = "Arthur Siegel <ajsiegel@optonline.com>"
0483__date__ = "Date: 2006-02-02"
0484__revision__ = "Revision: a1"
0485__url__ = "http://source.net/projects/pygeo"
0486__license__ ="GPL <http://www.opensource.org/licenses/gpl-license.php>"