Package topo :: Package coordmapper :: Module basic
[hide private]
[frames] | no frames]

Source Code for Module topo.coordmapper.basic

  1  """ 
  2  Functions for mapping from one 2D Cartesian coordinate system to another. 
  3   
  4  These are useful for defining projections with nonlinear magnifications, 
  5  reductions, or other transformations of the input coordinate system. 
  6   
  7  Generally, these function objects should work for an arbitrary (x,y) pair, 
  8  returning new, remapped (x,y), although some classes of mappers may 
  9  define their range and domain, with undesirable or undefined behavior 
 10  outside those regions.  It is best if such objects raise a suitable 
 11  exception in those circumstances. 
 12   
 13  $Id: basic.py 11296 2010-07-27 14:39:42Z ceball $ 
 14  """ 
 15  __version__='$Revision: 11296 $' 
 16   
 17  from math import atan,pi,atan2 
 18   
 19  from numpy import exp,log,sqrt,sin,cos,ones,dot 
 20  from numpy.matlib import matrix 
 21   
 22  import param 
 23   
 24  from topo.base.functionfamily import CoordinateMapperFn, IdentityMF 
 25  from topo.misc.util import signabs 
 26  from topo import numbergen 
 27   
 28   
 29   
 30   
31 -class ConstantMapper(CoordinateMapperFn):
32 """ 33 Map all values to the same constant, pre-specified coordinates. 34 """ 35 36 x_cons = param.Number(default=0.0,doc=""" 37 Constant x value returned by the mapping.""") 38 39 y_cons = param.Number(default=0.0, doc=""" 40 Constant y value returned by the mapping.""") 41
42 - def __call__(self, x, y):
43 # Ignores all (x,y), always returning (x_cons,y_cons) 44 return self.x_cons, self.y_cons
45 46
47 -class Pipeline(CoordinateMapperFn):
48 """ 49 Applies a sequence of coordmappers, left to right. 50 """ 51 52 mappers=param.List(default=[], 53 doc="The sequence of mappers to apply.") 54
55 - def __call__(self,x,y):
56 return reduce( lambda args,f: apply(f,args), 57 [(x,y)] + self.mappers )
58 59
60 -class Jitter(CoordinateMapperFn):
61 """ 62 Additively modifies calculated x,y coordinates, e.g. with random noise. 63 """ 64 65 scale = param.Parameter(default=0.0,doc="Amount of jitter.") 66 67 gen = param.Parameter(default=numbergen.UniformRandom(),doc= 68 "Number generator to use, typically a random distribution.") 69
70 - def __call__(self,x,y):
71 return x+(self.gen()-0.5)*self.scale,y+(self.gen()-0.5)*self.scale
72 73
74 -class NormalJitter(CoordinateMapperFn):
75 """ 76 Additively modifies calculated x,y coordinates with Gaussian random noise. 77 """ 78 79 gen = param.Parameter(default=numbergen.NormalRandom(),doc= 80 "Number generator to use, typically a Gaussian distribution.") 81
82 - def __call__(self,x,y):
83 return x+self.gen(),y+self.gen()
84 85
86 -class Grid(CoordinateMapperFn):
87 """ 88 Divides the 2D area into a grid, where all points within each grid 89 element map to the same location. 90 """ 91 92 xdensity = param.Number(default=1, bounds=(0,None), doc=""" 93 Number of columns per 1.0 input sheet distance horizontally.""") 94 95 ydensity = param.Number(default=1, bounds=(0,None), doc=""" 96 Number of rows per 1.0 input sheet distance vertically.""") 97
98 - def __call__(self,x,y):
99 xd=self.xdensity 100 yd=self.ydensity 101 102 xquant=(1.0/xd)*(int(xd*(x+0.5))-(0.5*(xd-1))) 103 yquant=(1.0/yd)*(int(yd*(y+0.5))-(0.5*(yd-1))) 104 105 return xquant,yquant
106 107
108 -class Polar2Cartesian(CoordinateMapperFn):
109 """ 110 Map from polar (radius,angle) to Cartesian (x,y) coordinates. 111 """ 112 113 degrees=param.Boolean(default=True, 114 doc="Indicates whether the input angle is in degrees or radians.") 115
116 - def __call__(self, r, theta):
117 118 if self.degrees: 119 theta = theta * pi/180 120 121 return r*cos(theta), r*sin(theta)
122 123
124 -class Cartesian2Polar(CoordinateMapperFn):
125 """ 126 Maps from Cartesian (x,y) to polar (radius,angle). 127 """ 128 129 degrees = param.Boolean(default=True, 130 doc="Indicates whether the output angle is in degrees or radians.") 131 132 negative_radii = param.Boolean(default = False, 133 doc="""If true, coordinates with negative x values will be 134 given negative radii, and angles between -90 and 90 135 degrees. (useful for mapping to saccade amplitude/direction space)""") 136
137 - def __call__(self, x, y):
138 139 if self.negative_radii: 140 xsgn,xabs = signabs(x) 141 radius = xsgn * sqrt(x*x+y*y) 142 angle = atan2(y,xabs) 143 else: 144 radius = sqrt(x*x+y*y) 145 angle = atan2(y,x) 146 147 if self.degrees: 148 angle *= 180/pi 149 150 return radius,angle
151 152 153 154
155 -class AffineTransform(CoordinateMapperFn):
156 """ 157 Remaps the input with an affine transform. 158 159 This mapper allows the specification of an arbitrary combination 160 of translation, rotation and scaling via a transform 161 matrix. Single translations, etc, can be specified more simply 162 with the subclasses Translate2d, Rotate2d, and Scale2d. 163 """ 164 165 matrix = param.Parameter(default=ones((3,3)),doc=""" 166 The affine transformation matrix. The functions 167 Translate2dMat, Rotate2dMat, and Scale2dMat generate affine 168 transform matrices that can be multiplied together to create 169 combination transforms. E.g. the matrix:: 170 171 Translate2dMat(3,0)*Rotate2d(pi/2) 172 173 will shift points to the right by 3 units and rotate them around 174 the origin by 90 degrees.""") 175
176 - def __init__(self, **kw):
177 super(AffineTransform,self).__init__(**kw) 178 179 # This buffer prevents having to allocate memory for each point. 180 self._op_buf = matrix([[0.0], 181 [0.0], 182 [1.0]])
183
184 - def __call__(self, x, y):
185 186 ## JPHACKALERT: If the coordmapper interface took a matrix of 187 ## x/y column vectors, instead of x and y separately, affine 188 ## transforms could be applied to all the points in a single 189 ## matrix operation. This would probably require revision of 190 ## some of the other coordmapper functions, but it might allow 191 ## some optimization in, e.g., CFProjection intialization by 192 ## allowing all the CF positions to be computed at once. 193 194 ## JPALERT: This is the easy way, but it allocates a matrix 195 ## for the result. It might be faster to unroll the 196 ## computation. 197 198 self._op_buf[0] = x 199 self._op_buf[1] = y 200 result = dot(self.matrix,self._op_buf) 201 202 return result[0,0],result[1,0]
203 204
205 -def Translate2dMat(xoff,yoff):
206 """ 207 Return an affine transformation matrix that translates points by 208 the offset (xoff,yoff). 209 """ 210 return matrix([[1, 0, xoff], 211 [0, 1, yoff], 212 [0, 0, 1 ]])
213 214
215 -def Rotate2dMat(t):
216 """ 217 Return an affine transformation matrix that rotates the points 218 around the origin by t radians. 219 """ 220 return matrix([[cos(t), -sin(t), 0], 221 [sin(t), cos(t), 0], 222 [ 0 , 0 , 1]])
223 224
225 -def Scale2dMat(sx,sy):
226 """ 227 Return an affine translation matrix that scales the points 228 toward/away from the origin by a factor of sx on x-axis and sy on 229 the y-axis. 230 """ 231 return matrix([[sx, 0, 0], 232 [ 0, sy, 0], 233 [ 0, 0, 1]])
234 235
236 -class Translate2d(AffineTransform):
237 """ 238 Translate the input by xoff,yoff. 239 """ 240 xoff = param.Number(default=0.0) 241 yoff = param.Number(default=0.0) 242
243 - def __init__(self,**kw):
244 super(Translate2d,self).__init__(**kw) 245 self.matrix = Translate2dMat(self.xoff,self.yoff)
246 247
248 -class Rotate2d(AffineTransform):
249 """ 250 Rotate the input around the origin by an angle in radians. 251 """ 252 angle = param.Number(default=0.0) 253
254 - def __init__(self,**kw):
255 super(Rotate2d,self).__init__(**kw) 256 self.matrix = Rotate2dMat(self.angle)
257 258
259 -class Scale2d(AffineTransform):
260 """ 261 Scale the input along the x and y axes by sx and sy, respectively. 262 """ 263 sx = param.Number(default=1.0) 264 sy = param.Number(default=1.0) 265
266 - def __init__(self,**kw):
267 super(Scale2d,self).__init__(**kw) 268 self.matrix = Scale2dMat(self.sx,self.sy)
269 270
271 -class SingleDimensionMapper(CoordinateMapperFn):
272 """ 273 Coordinate Mapper that uses an origin-centered 1-D mapping function. 274 275 An abstract mapping function for coordinate mappers that remap 276 based on the radius, x, or y individually. Subclasses should override 277 _map_fn(self,z). 278 """ 279 __abstract = True 280 281 in_range = param.Number(default=0.5*sqrt(2),bounds=(0,None),doc=""" 282 The maximum range of the mapping input.""") 283 out_range = param.Number(default=0.5*sqrt(2),bounds=(0,None), doc=""" 284 The maximum range of the output.""") 285 remap_dimension = param.ObjectSelector(default='radius', 286 objects=['radius','x','y','xy'],doc=""" 287 The dimension to remap. ('xy' remaps x and y independently.)""") 288 289
290 - def __call__(self,x,y):
291 292 if self.remap_dimension == 'radius': 293 r = sqrt(x**2 + y**2) 294 a = atan2(x,y) 295 new_r = self._map_fn(r) 296 xout = new_r * sin(a) 297 yout = new_r * cos(a) 298 else: 299 if 'x' in self.remap_dimension: 300 xout = self._map_fn(x) 301 else: 302 xout = x 303 304 if 'y' in self.remap_dimension: 305 yout = self._map_fn(y) 306 else: 307 yout = y 308 309 return xout,yout
310
311 - def _map_fn(self,z):
312 raise NotImplementedError
313 314
315 -class MagnifyingMapper(SingleDimensionMapper):
316 """ 317 Exponential (magnifying) mapping function. 318 319 Provides a mapping that magnifies the center of the activity image. 320 Parameter k indicates amount of magnification, where 0 means no 321 magnification. 322 """ 323 324 k = param.Number(default=1.0,bounds=(0,None)) 325
326 - def _map_fn(self,z):
327 k = self.k 328 if k == 0: 329 return z 330 else: 331 sgn,z = signabs(z) 332 return sgn * self.out_range * (exp(z/self.in_range*k)-1)/(exp(k)-1)
333 334
335 -class ReducingMapper(SingleDimensionMapper):
336 """ 337 Provides a mapping that reduces the center of the activity. 338 339 Roughly the inverse of MagnifyingMapper. k indicates amount of reduction. 340 """ 341 k = param.Number(default=1.0,bounds=(0,None)) 342
343 - def _map_fn(self,z):
344 k = self.k 345 sgn,z = signabs(z) 346 return sgn * self.out_range * log(z/self.in_range*k+1)/log(k+1)
347 348
349 -class OttesSCMapper(CoordinateMapperFn):
350 """ 351 Abstract class for superior colliculus mappings. 352 353 Subclasses of this class implement afferent and efferent mappings 354 from Ottes et al. (1986) Vision Research 26:857-873. 355 356 Default constant values are from Table 1, ibid. 357 """ 358 __abstract = True 359 360 361 A = param.Number(default=5.3, doc=""" 362 Shape parameter A, in degrees""") 363 Bu = param.Number(default=1.8, doc=""" 364 Rostral-caudal scale parameter, in mm""") 365 Bv = param.Number(default=1.8, doc=""" 366 Orthogonal (medial-lateral?) scale paraemter in mm/deg""") 367 368 mm_scale = param.Number(default=8.0,doc=""" 369 Scale factor to convert constants Bu and Bv from mm to sheet 370 units. Expressed in mm/unit. """) 371 372 amplitude_scale = param.Number(default=1,doc=""" 373 Scale factor for saccade command amplitude, expressed in 374 degrees per unit of sheet. Indicates how large a 375 saccade is represented by the x-component of the command 376 input.""") 377 378 direction_scale = param.Number(default=1,doc=""" 379 Scale factor for saccade command direction, expressed in 380 degrees per unit of sheet. Indicates what direction of saccade 381 is represented by the y-component of the command input.""") 382 383
384 - def __call__(self,x,y):
385 raise NotImplementedError
386 387
388 -class OttesSCMotorMapper(OttesSCMapper):
389 """ 390 Efferent superior colliculus mapping. 391 392 Provides the output/motor mapping from SC as defined by Ottes et al. 393 (see superclass docs for reference) 394 395 The mapping allows the creation of a single sheet representing 396 both colliculi, one in the x-positive hemisheet and the other in 397 the x-negative hemisheet. 398 """
399 - def __call__(self,x,y):
400 401 A = self.A 402 Bu = self.Bu / self.mm_scale 403 Bv = self.Bv / self.mm_scale 404 405 R = x * self.amplitude_scale 406 phi = y * self.direction_scale 407 408 Rsign,R = signabs(R) 409 410 u,v = ottes_mapping(R,phi,A,Bu,Bv) 411 return Rsign*u,v
412 413 414
415 -class OttesSCSenseMapper(OttesSCMapper):
416 """ 417 Afferent superior colliculus mapping. 418 419 Provides the input/retinal mapping to SC as defined by Ottes et al. 420 (see superclass docs for reference) 421 422 The mapping allows the creation of a single sheet representing 423 both colliculi, one in the x-positive hemisheet and the other in 424 the x-negative hemisheet. 425 426 [NOTE: see warning in docs for ottes_inverse_mapping().] 427 """ 428
429 - def __call__(self,x,y):
430 431 A = self.A 432 Bu = self.Bu 433 Bv = self.Bv 434 435 u = x * self.mm_scale 436 v = y * self.mm_scale 437 438 usgn,u = signabs(u) 439 440 R,phi = ottes_inverse_mapping(u,v,A,Bu,Bv) 441 442 return (usgn*R/self.amplitude_scale, 443 phi/self.direction_scale)
444 445
446 -def ottes_mapping(R,phi,A,Bu,Bv):
447 """ 448 The efferent mapping function from Ottes et al. (1986) 449 Vision Research 26:857-873. 450 451 Takes saccade with amplitude R (in degrees) and direction 452 phi (in degrees), and returns a location u,v on the colliculus 453 in mm, where the u axis is rostral/caudal, and the v axis is 454 medial/lateral. 455 """ 456 457 phi *= pi/180 458 u = Bu * (log(sqrt(R**2 + A**2 + 2*A*R*cos(phi))) - log(A)) 459 v = Bv * atan((R*sin(phi))/(R*cos(phi)+A)) 460 return u,v
461 462
463 -def ottes_inverse_mapping(u,v,A,Bu,Bv):
464 """ 465 Approximate inverse of ottes_mapping(), using the inverse function 466 provided in the appendix of Ottes et al. (1986) Vision Research 26:857-873 467 468 Takes takes a location u,v on the colliculus in mm and maps 469 to a retinal eccentricity R and direction phi, in degrees. 470 Inverse is approximate, with increasing error as positions near the 471 edges of the collicular sheet. (I.e. with high absolute v value). 472 """ 473 474 rads = pi/180 475 R = A * sqrt(exp(2*u/Bu) - 2*exp(u/Bu)*cos(rads*v/Bv) + 1) 476 #phi = atan( (exp(u/Bu)*sin(rads*v/Bv)) / (exp(u/Bu)*cos(rads*v/Bv) -1) ) 477 phi = atan2( (exp(u/Bu)*sin(rads*v/Bv)), (exp(u/Bu)*cos(rads*v/Bv) -1) ) * 180/pi 478 479 # JPALERT: Don't know why we have to multiply by 180/pi twice, but the answers 480 # are way off without it. Is the bug in my code, or in the original formula? 481 return R,phi*180/pi
482 483 484 # JPALERT: Temporary testing function. Will disappear eventually.
485 -def test_ottes_inverse():
486 487 A,Bu,Bv = 5.3,1.8,1.8 488 489 print '%10s %10s | %10s %10s | %10s %10s | %s' \ 490 % ('R in','phi in','R out','phi out','R err','phi err','phi ratio') 491 for r in range(10,60,10): 492 for phi in range(-60,60,10): 493 if phi != 0: 494 u,v = ottes_mapping(r,phi,A,Bu,Bv) 495 r2,phi2 = ottes_inverse_mapping(u,v,A,Bu,Bv) 496 print '%10.2f %10.2f | %10.2f %10.2f | %10.2f %10.2f | %.2f' \ 497 % (r,phi,r2,phi2,r2-r,phi2-phi,phi/phi2)
498 499 500 __all__ = list(set([k for k,v in locals().items() 501 if isinstance(v,type) \ 502 and issubclass(v,CoordinateMapperFn)])) 503