1 """
2 * Copyright (C) 2004, 2005 Cristina V. Torres
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
17 * MA 02111-1307 USA
18 """
19 __author__ = 'Cristina Torres <cristina.torres@ligo.org>'
20 __date__ = '$Date$'
21 __version__ = ''
22
23
24 import os
25 from numpy import ndarray,dtype,NaN,Inf,zeros,mean,sin,\
26 pi,all,any,array,arange,inner,cross,isnan,arcsin,\
27 mod
28 from pylal.metaarray import TimeSeries,TimeSeriesList
29 import copy
30 import sys
31 from StringIO import StringIO
32 from commands import getstatusoutput
33 from scipy import interp
34 from followup_utils import connectToFrameData
35 from pylal import metaarray
36 import warnings
37
38 try:
39 import sqlite3 as sqlite
40 except ImportError:
41 import sqlite
42 if os.getenv("DISPLAY") == None:
43
44 try:
45 import matplotlib
46 matplotlib.use("Agg")
47 import pylab
48 except Exception, errorInfo:
49 disableGraphics=True
50 sys.stderr.write("Error trying to import NON-INTERACTIVE pylab!\n")
51 sys.stderr.write("Exception Instance :%s\n"%(str(type(errorInfo))))
52 sys.stderr.write("Exception Args :%s\n"%(str(errorInfo.args)))
53 sys.stderr.write("Pylab functionality unavailable!\n")
54 else:
55
56 try:
57 import pylab
58 except Exception, errorInfo:
59 disableGraphics=True
60 sys.stderr.write("Error trying to import INTERACTIVE pylab!\n")
61 sys.stderr.write("Exception Instance :%s\n"%(str(type(errorInfo))))
62 sys.stderr.write("Exception Args :%s\n"%(str(errorInfo.args)))
63 sys.stderr.write("Pylab functionality unavailable!\n")
64
65 """
66 This file will hold all the methods required to run the c^2t(Cristina
67 Cesar Triangulation) DetChar tool.
68 """
69
70
71
72
73
75 """
76 Input two matrices Nx3 in size.
77 Calculates the inner product row by row.
78 """
79 if( A.shape != B.shape) and A.ndim > 2:
80 os.abort()
81 if A.ndim == 2:
82 return array([inner(A[i],B[i]) for i in arange(A.shape[0])])
83 else:
84 return array(inner(A,B))
85
86
87
89 """
90 Input two matrices Nx3 in size.
91 Calculates the cross product row by row.
92 """
93 if( A.shape != B.shape) and A.ndim > 2:
94 os.abort()
95 if A.ndim == 2:
96 return array([cross(A[i],B[i]) for i in arange(A.shape[0])])
97 else:
98 return array(cross(A,B))
99
100
101
102
104 """
105 Expects a 3xM matrices of 3D vectors for the plane norms
106 also in addition to that we take vectors from the "origin"
107 (0,0,0) (single 3x1 matrix) to be the coordinate input. The
108 output is a tuple of 3xM arrays parameterized by their index,
109 intersection vector and point on that line. Output is
110 (linePoint,LineVector) ((3xM),(3xM)). The input wiggleA is the
111 wiggle angle in radian to consider two vectors parrallel.
112 """
113
114
115 for j in arange(len(normA)):
116 normA[j]=normA[j]/CSTdot(normA[j],normA[j])
117 normB[j]=normB[j]/CSTdot(normB[j],normB[j])
118
119 lineVec=zeros(normA.shape)
120 linePoi=zeros(normA.shape)
121 for i in arange(len(normA)):
122
123 if arcsin(CSTdot(normA[i],normB[i])) >= (1-wiggleA):
124
125 lineVec[i]=NaN
126 linePoi[i]=NaN
127 else:
128
129 tmpVec=CSTcross(normA[i],normB[i])
130 lineVec[i]=tmpVec/CSTdot(tmpVec,tmpVec)
131 print "Intersecting line vector(unit) %s"%(lineVec)
132
133 if isnan(lineVec[i]).any() or all(lineVec[i]==0.0):
134 print "Vector of plane intersection line is ZERO Vector!"
135 lineVec[i]=NaN
136 linePoi[i]=NaN
137 else:
138
139 zeroIndex=lineVec.argmax()
140
141 (iX,iY,iZ)=mod(range(zeroIndex,zeroIndex+3),3).tolist()
142 print "Vector elements:%s,%s,%s"%(iX,iY,iZ)
143
144 linePoi[i][iZ]=0
145 linePoi[i][iX]=(CSTdot(normB[i],coordB)-(CSTdot(normA[i],coordA)*normB[i][iY])/normA[i][iY])/(normB[i][iX]-(normB[i][iY]*normA[i][iX])/normA[i][iY])
146 linePoi[i][iY]=(CSTdot(normA[i],coordA)-normA[i][iX]*linePoi[i][iX])/normA[i][iY]
147 return (lineVec,linePoi)
148
149
150
151
153 """
154 Expects inputs of 3xM matrices representing vectors along a line
155 and a vector 3xM of matching points on lines for those vectors.
156 Output is in a 3xM collection of intersecting points.
157 Examples:
158 Input
159 ([[vX,vY,vZ]],[[lX,lY,lZ]])
160
161 Output
162 [[i,j,k]]
163 A set of lines is considered parallel if angle(rads) between
164 then is less than wiggleA, in this case the output coordinate
165 is (NaN,NaN,NaN).
166 """
167
168
169 for j in arange(len(vecA)):
170 vecA[j]=vecA[j]/CSTdot(vecA[j],vecA[j])
171 vecB[j]=vecB[j]/CSTdot(vecB[j],vecB[j])
172
173 triPoi=zeros(vecA.shape)
174 for i in arange(len(vecA)):
175
176 if arcsin(CSTdot(vecA[i],vecB[i]))<=wiggleA:
177
178 triPoi[i]=NaN
179 else:
180 (iX,iY,iZ)=(0,1,2)
181
182
183 gammaConst=(\
184 coordA[i][iY]*vecA[i][iZ] - vecA[i][iZ]*coordB[i][iY] +\
185 vecA[i][iY]*(coordB[i][iZ]-coordA[i][iZ])\
186 )/\
187 (vecA[i][iY]*vecB[i][iZ]+vecA[i][iZ]*vecB[i][iY])
188 triPoi=coordB+(gammaConst*vecB)
189 return triPoi
190
192 """
193 This class is just method globber that works expects as input
194 TimeSeries objects defined by pylal.frutils. This class only
195 works in a cartesian basis! The distance should be measure in units
196 of meters.
197 """
198 - def __init__(self,\
199 name="Unknown",\
200 xAxis=None,\
201 yAxis=None,\
202 zAxis=None,\
203 coordinate=(None,None,None),\
204 swapSign=(False,False,False),\
205 tType="normal"\
206 ):
207 """
208 This object is initialized with 3 data streams (1xM) which should be the
209 intesities of each spatial coordinate seen by the sensor in
210 question. In addition to this information, the cartesian
211 coordinates relative to some origin should be specified. Inside
212 the class the data will be represented as collections of row
213 vectors in a 3xM matrix. The units
214 are not checked between different cSqrTSensor objects. The user
215 should take care to make sure this is sensible!
216 """
217 self.name=name.strip().upper()
218 self.__tType__=tType.strip().lower()
219 if self.__tType__ != "normal":
220 os.abort()
221
222
223 self.dtype='float64'
224
225
226 wiggleAngle=3.0
227 wiggleDistance=25.0
228 self.__wiggleA__=None
229 self.__wiggleD__=None
230 self.setParallelMismatchAngle(wiggleAngle)
231 self.setHypersphereMismatch(wiggleDistance)
232
233
234
235 if isinstance(xAxis,TimeSeries):
236 self.metaX=xAxis.metadata
237 else:
238 self.metaX=None
239 if isinstance(yAxis,TimeSeries):
240 self.metaY=yAxis.metadata
241 else:
242 self.metaY=None
243 if isinstance(zAxis,TimeSeries):
244 self.metaZ=zAxis.metadata
245 else:
246 self.metaZ=None
247
248
249
250 if all(self.metaX,self.metaY,self.metaZ):
251 if not (
252 xMetaDict["segments"].intersects_all(yMetaDict["segments"]) \
253 and \
254 yMetaDict["segments"].intersects_all(zMetaDict["segments"]) \
255 and \
256 zMetaDict["segments"].intersects_all(xMetaDict["segments"]) \
257 ):
258 raise Exception, \
259 "Inputs for X,Y,Z data don't cover identical segments!"
260
261
262
263 if not(\
264 xMetaDict["dt"] == \
265 yMetaDict["dt"] == \
266 zMetaDict["dt"] \
267 ):
268 raise Exception, \
269 "Input data for sampling rates on data is not consistent!"
270 else:
271 sys.stdout.write("Inputs are not TimeSeries ojects can not verify data traits.\n")
272
273
274
275
276 if isinstance(coordinate,ndarray):
277 self.coordinate=coordinate
278 else:
279 self.coordinate=array(coordinate,dtype=self.dtype)
280
281
282
283 self.data=zeros((len(xAxis),3),self.dtype)
284 self.data[:,0]=xAxis
285 self.data[:,1]=yAxis
286 self.data[:,2]=zAxis
287
288
289
290
294 """
295 Enter an angle in degrees to define a tolerance on the
296 parallelness of two lines for determining potential
297 intersection of lines and planes.
298 """
299 self.__wiggleA__=sin(wiggleAngle/(2.0*pi))
300
304 """
305 Enter an angle in degrees to define a tolerance on the
306 parallelness of two lines for determining potential
307 intersection of lines and planes.
308 """
309 self.__wiggleD__=wiggleDistance
310
311 - def triangulate(self,\
312 otherSensor=[None],\
313 ):
314 """
315 Wrapper to pick the correct triangulation machinery.
316 """
317 if not isinstance(otherSensor,list):
318 otherSensor=[otherSensor]
319 if not all([isinstance(x,cSqrTSensor) for x in otherSensor]):
320 raise Exception, "Expected list of cSqrTSensor objects, got list of other objects!"
321 if self.__tType__ == "normal":
322 return self.__normalTriangulate__(otherSensor)
323 else:
324 os.abort()
325
329 """
330 Assuming the input is normals to 3D planes in 3 space we will
331 intercept the planes resulting in either points in 3D space or 2D
332 lines whichever is the most physically interesting.
333 """
334
335
336 if self.name in [x.name for x in otherSensor]:
337 raise Exception, "Namespace collision between sensors to triangulate!"
338
339
340
341 planeIntersections=dict()
342 pKeyMask="%s_%s"
343 if len(otherSensor) >= 1:
344
345 sensorList=list()
346 sensorList.append(self)
347 sensorList.extend(otherSensor)
348
349 for aa,A in enumerate(sensorList):
350 for bb,B in enumerate(sensorList):
351 myKey=pKeyMask%(A.name,B.name)
352 symKey=pKeyMask%(B.name,A.name)
353 if (aa!=bb) and \
354 not (myKey in planeIntersections.keys()) and \
355 not (symKey in planeIntersections.keys()):
356 planeIntersections[pKeyMask%(A.name,B.name)]=intersectPlanes(A.data,\
357 B.data,\
358 A.coordinate,\
359 B.coordinate,
360 wiggleA=self.__wiggleA__)
361 else:
362 raise Exception, "Need at least two normals to intersect planes."
363
364 if len(planeIntersections.keys()) == 1:
365 myKey=planeIntersections.keys()
366 return planeIntersections[myKey[0]]
367
368 lineIntersections=dict()
369 for (nameA,lineA) in planeIntersections.iteritems():
370 for (nameB,lineB) in planeIntersections.iteritems():
371 lineIntersections[pKeyMask%(nameA,nameB)]=intersectLines(lineA[0],\
372 lineB[0],\
373 lineA[1],\
374 lineB[1],\
375 wiggleA=self.__wiggleA__)
376
377
378
379 resultPoint=zeros(lineIntersections[0],dtype=self.dtype)
380 resultVectors=zeros(lineIntersections[0],dtype=self.dtype)
381
382 for i in arange(len(lineIntersections[0])):
383 tmpCentroid=zeros(lineIntersections[0],dtype=self.dtype)
384 for j in arange(len(lineIntersections)):
385 tmpCentroid+=lineIntersections[j]
386 tmpCentroid=tmpCentroid/float(len(lineIntersections))
387 resultPoint[i]=tmpCentroid
388 if max([sqrt(CSTdot(x-tmpCentroid,x-tmpCentroid)) \
389 for x in lineIntersections.iteritems]) > \
390 self.__wiggleD__:
391 resultPoint[i]=NaN
392 else:
393 resultPoint[i]=tmpCentroid
394
395
396
397 return (resultVectors,resultPoint)
398