1 |
# $Id:$ |
2 |
""" |
3 |
transformations |
4 |
|
5 |
@var __author__: name of author |
6 |
@var __copyright__: copyrights |
7 |
@var __license__: licence agreement |
8 |
@var __url__: url entry point on documentation |
9 |
@var __version__: version |
10 |
@var __date__: date of the version |
11 |
@var DEG: unit of degree |
12 |
@var RAD: unit of radiant |
13 |
""" |
14 |
|
15 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
16 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
17 |
http://www.access.edu.au |
18 |
Primary Business: Queensland, Australia""" |
19 |
__license__="""Licensed under the Open Software License version 3.0 |
20 |
http://www.opensource.org/licenses/osl-3.0.php""" |
21 |
__url__="http://www.iservo.edu.au/esys/escript" |
22 |
__version__="$Revision:$" |
23 |
__date__="$Date:$" |
24 |
|
25 |
import numarray |
26 |
import math |
27 |
|
28 |
_TYPE=numarray.Float64 |
29 |
DEG=math.pi/180. |
30 |
RAD=1. |
31 |
class Transformation(object): |
32 |
""" |
33 |
general class to define an affine transformation x->Ax+b |
34 |
""" |
35 |
def __init__(self): |
36 |
""" |
37 |
create a linear transformation |
38 |
""" |
39 |
pass |
40 |
def __call__(self,x=numarray.zeros((3,))): |
41 |
""" |
42 |
applies transformation to x |
43 |
""" |
44 |
raise NotImplementeError() |
45 |
class Translation(Transformation): |
46 |
""" |
47 |
defines a translation x->x+b |
48 |
""" |
49 |
def __init__(self,b=numarray.zeros((3,),type=_TYPE)): |
50 |
""" |
51 |
create linear transformation x->x+b |
52 |
""" |
53 |
super(Translation, self).__init__() |
54 |
self.__b=numarray.array(b,_TYPE) |
55 |
def __call__(self,x=numarray.zeros((3,))): |
56 |
""" |
57 |
applies translation to x |
58 |
""" |
59 |
return numarray.array(x,_TYPE)+self.__b |
60 |
class Rotatation(Transformation): |
61 |
""" |
62 |
defines a rotation |
63 |
""" |
64 |
def __init__(self,point0=numarray.zeros((3,),type=_TYPE),point1=numarray.ones((3,),type=_TYPE),angle=0.*RAD): |
65 |
""" |
66 |
creates a rotation using to points to define the axis and a rotation angle |
67 |
""" |
68 |
self.__point0=numarray.array(point0,type=_TYPE) |
69 |
self.__point1=numarray.array(point1,type=_TYPE) |
70 |
self.__axis=self.__point1-self.__point0 |
71 |
lax=numarray.dot(self.__axis,self.__axis) |
72 |
if not lax>0: |
73 |
raise ValueError("points must be distinct.") |
74 |
self.__axis/=math.sqrt(lax) |
75 |
self.__angle=float(angle) |
76 |
def __call__(self,x=numarray.zeros((3,))): |
77 |
""" |
78 |
applies rotatation to x |
79 |
""" |
80 |
x=numarray.array(x,_TYPE) |
81 |
z=x-self.__point0 |
82 |
z0=numarray.dot(z,self.__axis) |
83 |
z_per=z-z0*self.__axis |
84 |
lz_per=numarray.dot(z_per,z_per) |
85 |
if lz_per>0: |
86 |
axis1=z_per/math.sqrt(lz_per) |
87 |
axis2=_cross(axis1,self.__axis) |
88 |
lax2=numarray.dot(axis2,axis2) |
89 |
if lax2>0: |
90 |
axis2/=math.sqrt(lax2) |
91 |
return z0*self.__axis+math.sqrt(lz_per)*(math.cos(self.__angle)*axis1-math.sin(self.__angle)*axis2)+self.__point0 |
92 |
else: |
93 |
return x |
94 |
else: |
95 |
return x |
96 |
def _cross(x, y): |
97 |
""" |
98 |
Returns the cross product of x and y |
99 |
""" |
100 |
return numarray.array([x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0]], _TYPE) |
101 |
|
102 |
class Dilation(Transformation): |
103 |
""" |
104 |
defines a dilation |
105 |
""" |
106 |
def __init__(self,factor=1.,center=numarray.zeros((3,),type=_TYPE)): |
107 |
""" |
108 |
creates a dilation with a center an a given expansion/contraction factor |
109 |
""" |
110 |
if not factor>0: |
111 |
raise ValueError("factor must be positive.") |
112 |
self.__factor=factor |
113 |
self.__center=numarray.array(center,type=_TYPE) |
114 |
def __call__(self,x=numarray.zeros((3,))): |
115 |
""" |
116 |
applies dilation to x |
117 |
""" |
118 |
x=numarray.array(x,_TYPE) |
119 |
return self.__factor*(x-self.__center)+self.__center |
120 |
|
121 |
class Reflection(Transformation): |
122 |
""" |
123 |
defines a reflection on a plain |
124 |
""" |
125 |
def __init__(self,normal=numarray.ones((3,),type=_TYPE),offset=0.): |
126 |
""" |
127 |
defines a reflection on a plain defined in normal form |
128 |
""" |
129 |
self.__normal=numarray.array(normal,type=_TYPE) |
130 |
ln=math.sqrt(numarray.dot(self.__normal,self.__normal)) |
131 |
if not ln>0.: |
132 |
raise ValueError("normal must have positive length.") |
133 |
self.__normal/=ln |
134 |
if isinstance(offset,float) or isinstance(offset,int): |
135 |
self.__offset=offset/ln |
136 |
else: |
137 |
self.__offset=numarray.dot(numarray.array(offset,type=_TYPE),self.__normal) |
138 |
def __call__(self,x=numarray.zeros((3,))): |
139 |
""" |
140 |
applies reflection to x |
141 |
""" |
142 |
x=numarray.array(x,_TYPE) |
143 |
return x - 2*(numarray.dot(x,self.__normal)-self.__offset)*self.__normal |