WikiStart: geometry.py

File geometry.py, 20.4 kB (added by jbosch, 1 year ago)

tutorial-enhanced ellipse class

Line 
1 import numpy
2 import weakref
3 from util import AttributeDict, BoundClass
4
5 # Set this to 1.0 to use radians instead of degrees.
6 ANGLE_UNITS = numpy.pi/180.0
7
8 class Ellipse(object):
9     """Represents an ellipse with one of the following parameterizations:
10      - center x,y + second moments (xx,yy,xy)
11      - center x,y + axes (a,b) + angle (theta)
12      - center x,y + complex ellipticity (e1,e2) or e + geometric mean radius (r)
13
14      The Ellipse class does not require that the units of the center match
15      those of the size, but other functions and classes that use Ellipse objects
16      may, so it might be a good idea to do that whenever possible..
17      """
18
19     def __init__(self,x=0.0,y=0.0,xx=1.0,yy=1.0,xy=0.0,moments=None,center=None):
20         # *TUTORIAL*
21         #
22         # I've found that it's often a good idea to make the __init__
23         # function accept all the arguments that make up the "true"
24         # attributes of the object.  You can add other optional
25         # arguments too, but often times that can be better handled by
26         # the 'named constructor' idiom - a static method that makes
27         # new objects in a special manner.
28         #
29         if center is not None:
30             x,y = center
31         self.x = float(x)
32         self.y = float(y)
33         if moments is not None:
34             xx,yy,xy = moments
35         self.xx = xx
36         self.yy = yy
37         self.xy = xy
38
39     @staticmethod
40     def from_moments(x=0.0,y=0.0,xx=None,yy=None,xy=None,moments=(1,1,0),center=None):
41         """Construct an Ellipse from centroid and second image moments."""
42         # *TUTORIAL*
43         #
44         # This is a static method, which means you can call it on the
45         # type ("Ellipse") or an instance of that type, and it won't
46         # pass a "self" object along.  It's basically just like a
47         # regular function, only it's stored inside a class.  The
48         # '@staticmethod' that makes this definition into a static
49         # method is called a "decorator".  There are a few other
50         # predefined decorators, and you can even make your own to do
51         # other interesting things to functions (and classes).
52         #
53         # This static method also happens to be a named constructor.
54         #
55         return Ellipse(x=x,y=y,xx=xx,yy=yy,xy=xy,moments=moments,center=center)
56
57     @staticmethod
58     def from_axes(x=0.0,y=0.0,a=1.0,b=1.0,theta=0.0,center=None):
59         """Construct an Ellipse from centroid, semi-major/minor axes, and position angle."""
60         self = Ellipse(x=x,y=y,center=center)
61         self.set_axes(a,b,theta)
62         return self
63
64     @staticmethod
65     def from_ellipticity(x=0.0,y=0.0,e1=0.0,e2=0.0,r=1.0,e=None,center=None):
66         """Construct an Ellipse from centroid, complex ellipticity, and radius."""
67         self = Ellipse(x=x,y=y,center=center)
68         self.set_ellipticity(e1,e2,r,e)
69         return self
70    
71     # *TUTORIAL*
72     #
73     # The next two function definitions and the following declaration
74     # are a "getter", a "setter", and a property.  This is a way to define
75     # calculated attributes, or possibly to put constraints on the values
76     # of regular attributes.  The property declaration takes those two
77     # functions and uses them to figure out what to do when a user
78     # tries to manipulate the "center" attribute of an ellipse.
79     #
80     # On Private Methods and Attributes:
81     #
82     # - The double underscores before the function names (and not
83     #   after them) indicate that __get_center and __set_center are
84     #   "private" to the Ellipse class: users shouldn't use them (only
85     #   the "center" attribute itself).
86     #
87     # - This makes it possible to later change how the center
88     #   attribute is implemented; maybe we later decide it's more
89     #   efficient to actually set "self.center=(x,y)" in __init__ - if
90     #   users depended on __set_center, their code that used Ellipse
91     #   objects would no longer work, but if they just always accessed
92     #   the center attribute, they wouldn't even notice the change.
93     #
94     # - It's not possible to make an attribute or method truly private
95     #   in Python, like it is in C++ or Java, but a single preceeding
96     #   underscore is a very standard convention that tells users
97     #   "don't use this unless you know what you are doing!".  The
98     #   double underscore is even stronger - it causes the name of the
99     #   attribute to be "mangled" outside the class.  You'll notice
100     #   that the Ellipse class, when viewed from outside, doesn't
101     #   appear to have a "__get_center" method.  But it does have a
102     #   "_Ellipse__get_center" method!
103     #
104     # - Name mangling can be confusing when you have a bunch of
105     #   related classes that need to access each other's private data.
106     #   In this case, it's much clearer to just use a single
107     #   preceeding underscore.  Avoid trying to outwit the name
108     #   mangling; even though it isn't hard to do.
109     #
110     # For all you experienced OOP programmers: properties aren't virtual.
111     #
112     # - Be aware: if you derive a new class based on Ellipse, and
113     #   override __get_center or __set_center (this is actually a
114     #   fairly common thing to do), the center property will still use
115     #   the old ones!
116     #
117     def __get_center(self): return (self.x,self.y)
118     def __set_center(self,center): self.x,self.y = center
119     center = property(__get_center,__set_center,doc="center point as a tuple")
120
121     def __get_moments(self): return (self.xx,self.yy,self.xy)
122     def __set_moments(self,moments): self.xx,self.yy,self.xy = moments
123     moments = property(__get_moments,__set_moments,doc="tuple of moments (xx,yy,xy)")
124
125     def set_axes(self,a=None,b=None,theta=None):
126         a_i,b_i,theta_i = self._compute_axes_from_moments(self.xx,self.yy,self.xy)
127         if a is None: a = a_i
128         if b is None: b = b_i
129         if theta is None: theta = theta_i
130         self.xx,self.yy,self.xy = self._compute_moments_from_axes(a,b,theta)
131
132     def set_ellipticity(self,e1=None,e2=None,r=None,e=None):
133         e1_i,e2_i,r_i = self._compute_ellipticity_from_moments(self.xx,self.yy,self.xy)
134         if e is not None:
135             e1 = e.real
136             e2 = e.imag
137         if e1**2 + e2**2 > 1.0:
138             raise ValueError("complex ellipticity must have magnitude <= 1.0")
139         if e1 is None: e1 = e1_i
140         if e2 is None: e2 = e2_i
141         if r is None: r = r_i
142         self.xx,self.yy,self.xy = self._compute_moments_from_ellipticity(e1,e2,r)
143
144     def __get_a(self): return self._compute_axes_from_moments(self.xx,self.yy,self.xy)[0]
145     def __set_a(self,a): self.set_axes(a=a)
146     a = property(__get_a,__set_a,doc="semi-major axis")
147
148     def __get_b(self): return self._compute_axes_from_moments(self.xx,self.yy,self.xy)[1]
149     def __set_b(self,b): self.set_axes(b=b)
150     b = property(__get_b,__set_b,doc="semi-minor axis")
151    
152     def __get_theta(self): return self._compute_axes_from_moments(self.xx,self.yy,self.xy)[2]
153     def __set_theta(self,theta): self.set_axes(theta=theta)
154     theta = property(__get_theta,__set_theta,doc="position angle of major axis, CCW from x-axis")
155
156     def __get_e1(self): return self._compute_ellipticity_from_moments(self.xx,self.yy,self.xy)[0]
157     e1 = property(__get_e1,doc="real component of complex ellipticity, |e|=(a^2-b^2)/(a^2+b^2)")
158
159     def __get_e2(self): return self._compute_ellipticity_from_moments(self.xx,self.yy,self.xy)[1]
160     e2 = property(__get_e2,doc="imaginary component of complex ellipticity, |e|=(a^2-b^2)/(a^2+b^2)")
161
162     def __get_e(self):
163         e1,e2,r = self._compute_ellipticity_from_moments(self.xx,self.yy,self.xy)
164         return complex(e1,e2)
165     def __set_e(self,e): self.set_ellipticity(e=e)
166     e = property(__get_e,__set_e,doc="complex ellipticity, |e|=(a^2-b^2)/(a^2+b^2)")
167    
168     def __get_r(self): return self._compute_ellipticity_from_moments(self.xx,self.yy,self.xy)[2]
169     def __set_r(self,r): self.set_ellipticity(r=r)
170     r = property(__get_r,__set_r,doc="geometric mean radius, r = (a*b)^(1/2) = (xx*yy-xy^2)^(1/4)")
171
172     def convolve(self,other):
173         """Convolve two elliptical Gaussians represented by Ellipse objects,
174         using the centroid of self and ignoring the centroid of other.
175         """
176         return Ellipse.from_moments(self.x,self.y,self.xx+other.xx,self.yy+other.yy,self.xy+other.xy)
177
178     def __str__(self):
179         # *TUTORIAL*
180         #
181         # If you haven't seen dictionary-based string formatting
182         # before, it's pretty cool.  Unfortunately, we often want to
183         # format the attributes of an object, not the values from a
184         # dictionary (and in Python 3.0, there's an easy way to do
185         # this).  But for now, this is where the AttributeDict class
186         # comes in - it takes an object, and acts like a dictionary
187         # that turns 'obj["name"]' requests into 'obj.name'.
188         #
189         # The source code for AttributeDict is in the util module.
190         #
191         return "Ellipse(x=%(x)s,y=%(y)s,xx=%(xx)s,yy=%(yy)s,xy=%(xy)s,\n" \
192                "        a=%(a)s,b=%(b)s,theta=%(theta)s,\n" \
193                "        e1=%(e1)s,e2=%(e2)s,r=%(r)s)" % AttributeDict(self)
194
195     def __repr__(self):
196         return "Ellipse(x=%(x)r,y=%(y)r,xx=%(xx)r,yy=%(yy)r,xy=%(xy)r)" % AttributeDict(self)
197
198     def __reduce__(self):
199         # *TUTORIAL*
200         #
201         # If you want your custom object to be picklable, there's
202         # nothing special to do; it happens automatically.  If you
203         # want to control exactly how that happens, or you have some
204         # special attributes that shouldn't be pickled (such as the
205         # only sometimes-present matplotlib attribute of Ellipse), you
206         # can define a __reduce__ method, which is expected to return
207         # a callable object and a tuple of arguments to pass to it -
208         # the new object will be reconstructed from the return value.
209         # 99% of the time that callable should be the class itself,
210         # which implicitly calls __init__.  And that's why I like
211         # having __init__ be closely tied to the actual attributes of
212         # an object.
213         #
214         return Ellipse,(self.x,self.y,self.xx,self.yy,self.xy)
215
216     class RadialFunctor(object):
217         """A callable object that returns the distance from the ellipse center to another point,
218         divided by the radius of the ellipse along the line between them.
219         """
220         # *TUTORIAL*
221         #
222         # RadialFunctor is a nested class; it is defined inside
223         # another class, Ellipse.
224         # Nested classed aren't a big deal in Python, and they
225         # don't have any special privileges.  It's only useful
226         # to nest classes purely for organizational purposes.
227         #
228         # On caveat with nested classes is that they can't be
229         # pickled (this is just a limitation of how pickle
230         # works), so if you want a class to be picklable,
231         # don't make it a nested class.
232
233         def __init__(self,ellipse):
234             self.ellipse = ellipse
235             self.r_inv = 1.0/(self.ellipse.xx*self.ellipse.yy-self.ellipse.xy**2)
236
237         def __call__(self,x,y):
238             # *TUTORIAL*
239             #
240             # An instance of RadialFunctor is a callable object,
241             # often called a "functor".  It may have data and
242             # other methods, but it most importantly has a
243             # __call__ method, which defines what happens
244             # when you do, for instance, "rf(5,3)", when
245             # rf is an instance of RadialFunctor.
246             #
247             # Objects-as-functions is an incredibly powerful
248             # idea, and it leads to a whole new programming
249             # paradigm ("functional programming").
250             #
251             dx = x-self.ellipse.x
252             dy = y-self.ellipse.y
253             return numpy.sqrt(self.r_inv*(self.ellipse.yy*dx**2
254                                           - 2.0*dx*dy*self.ellipse.xy
255                                           + self.ellipse.xx*dy**2))
256
257         def grid(self,x,y):
258             """Evaluate the functor on a grid; equivalent to self(*numpy.meshgrid(x,y))
259             (but more efficient)."""
260             dx = x-self.ellipse.x
261             dy = y-self.ellipse.y
262             return numpy.sqrt(self.r_inv*(self.ellipse.yy*dx[numpy.newaxis,:]**2
263                                           - 2.0*numpy.multiply.outer(dy,dx)*self.ellipse.xy
264                                           + self.ellipse.xx*dy[:,numpy.newaxis]**2))
265
266     # *TUTORIAL*
267     #
268     # This next declaration looks a lot like a property declaration, but
269     # it uses something I created myself in the util module, "BoundClass".
270     # "BoundClass" and "property" instances are 'descriptors', a special type
271     # of Python object that is added as a class attribute and controls how
272     # to get an attribute of the same name from an instance.  Back up and
273     # read that again if you need to.
274     #
275     # Here's how it works:
276     #
277     # - The person writing a class (like Ellipse), adds a descriptor object
278     #   (an instance of BoundClass, called "radial"), to the class.
279     #
280     # - You access the BoundClass itself when you do Ellipse.radial (you're
281     #   accessing a class attribute).
282     #
283     # - When a a user asks for e.radial, when "e" is an Ellipse instance,
284     #   magic happens, and Python calls the __get__ attribute of radial,
285     #   (and similarly for setting e.radial and __set__).  For BoundClass,
286     #   this just returns a new RadialFunctor object initialized by e.
287     #   For a property, this would access the getter and setter functions.
288     #
289     # - Most of the time you'll just use property.  But it's nice to know
290     #   how it's working under the hood, and someday you may want to write
291     #   your own custom descriptor.
292     #
293     # - Note that gives us now four ways of customizing attribute access
294     #   in Python:
295     #   * define "real" attributes in __init__ (or anywhere else, actually)
296     #   * define functions inside the class definition (this makes methods)
297     #   * write __getattr__ and __setattr__ methods
298     #   * use descriptors
299     #
300     radial = BoundClass(RadialFunctor)
301
302     class ParametricFunctor(object):
303         """A callable object that returns a tuple of x,y for the ellipse border, given
304         a parameter t that runs from 0 to 2pi.
305         """
306
307         def __init__(self,ellipse):
308             self.x = ellipse.x
309             self.y = ellipse.y
310             angle = ellipse.theta*ANGLE_UNITS
311             c = numpy.cos(angle)
312             s = numpy.sin(angle)
313             a = ellipse.a
314             b = ellipse.b
315             self.ux = a*c
316             self.uy = a*s
317             self.vx = -b*s
318             self.vy = b*c
319
320         def __call__(self,t):
321             c = numpy.cos(t)
322             s = numpy.sin(t)
323             return (self.x + self.ux*c + self.vx*s, self.y + self.uy*c + self.vy*s)
324
325     parametric = BoundClass(ParametricFunctor)
326
327     # *TUTORIAL*
328     #
329     # The matplotlib interface for Ellipse is another
330     # nested class, but in this case rather than accessing
331     # it with a BoundClass that creates a temporary object
332     # that is usually thrown away, calling "plot()" on an
333     # ellipse object creates a new attribute "matplotlib"
334     # that sticks around so you can update the plot with
335     # different plotting styles or with updated ellipse
336     # parameters.  You can delete it with "del e.matplotlib"
337     # or replace it by calling "plot()" again at any time.
338     #
339     # There's some fairly fancy matplotlib programming
340     # in here if you're curious, but it won't be covered
341     # it this tutorial.
342     #
343     class MatplotlibInterface(object):
344         """An interface for drawing the Ellipse using matplotlib.
345
346         This is typically initiated by calling Ellipse.plot(), which
347         adds the interface as the matplotlib attribute of the ellipse
348         object (this can be deleted later if desired).
349         """
350        
351         def __init__(self,ellipse,points=None,**kwds):
352             import matplotlib.patches
353             self.__ellipse = weakref.proxy(ellipse)
354             if points is None:
355                 self.patch = matplotlib.patches.Ellipse((self.__ellipse.x,self.__ellipse.y),
356                                                          self.__ellipse.a,self.__ellipse.b,
357                                                          self.__ellipse.theta*ANGLE_UNITS*180.0/numpy.pi,
358                                                          **kwds)
359             else:
360                 t = numpy.linspace(0,numpy.pi*2.0,points)
361                 xy = numpy.zeros((points,2),dtype=float)
362                 xy[:,0],xy[:,1] = self.__ellipse.parametric(t)
363                 self.patch = matplotlib.patches.Polygon(xy,**kwds)
364
365         def __getattr__(self,name):
366             return getattr(self.patch,name)
367
368         def update(self,show=True,rescale=True):
369             """Update the matplotlib representation to the
370             current ellipse parameters.
371             """
372             import matplotlib.patches
373             if isinstance(self.patch,matplotlib.patches.Polygon):
374                 points = self.patch.get_xy().shape[0]
375                 t = numpy.linspace(0,numpy.pi*2.0,points)
376                 xy = numpy.zeros((points,2),dtype=float)
377                 xy[:,0],xy[:,1] = self.__ellipse.parametric(t)
378                 self.patch.set_xy(xy)
379                 axes = self.patch.get_axes()
380             elif isinstance(self.patch,matplotlib.patches.Ellipse):
381                 new_patch =  matplotlib.patches.Ellipse((self.__ellipse.x,self.__ellipse.y),
382                                                          self.__ellipse.a,self.__ellipse.b,
383                                                          self.__ellipse.theta*ANGLE_UNITS*180.0/numpy.pi)
384                 new_patch.update_from(self.patch)
385                 axes = self.patch.get_axes()
386                 if axes is not None:
387                     self.patch.remove()
388                     axes.add_patch(new_patch)
389                 self.patch = new_patch
390             if axes is not None:
391                 if rescale: axes.autoscale_view()
392                 if show: axes.figure.canvas.draw()
393
394     def plot(self,axes=None,show=True,rescale=True,points=None,**kwds):
395         """Plot the Ellipse in matplotlib, adding a MatplotlibInterface
396         object as the 'matplotlib' attribute of the ellipse.
397         
398         Aside from those below, keyword arguments for the
399         matplotlib.patches.Patch constructor are also accepted
400         ('facecolor', 'linestyle', etc.)
401
402         Arguments:
403            axes -------- A matplotlib.axes.Axes object, or None to use
404                          matplotlib.pyplot.gca().
405            show -------- If True, update the figure automatically.  Set
406                          to False for batch processing.
407            rescale ----- If True, rescale the axes.
408            points ------ Number of points to use in approximating the
409                          Ellipse.  If None, use a (slower?) scale-free
410                          representation.
411         """
412         import matplotlib.pyplot
413         self.matplotlib = Ellipse.MatplotlibInterface(self,points,**kwds)
414         if axes is None:
415             axes = matplotlib.pyplot.gca()
416         axes.add_patch(self.matplotlib.patch)
417         if rescale: axes.autoscale_view()
418         if show: axes.figure.canvas.draw()
419         return self.matplotlib.patch           
420
421     @staticmethod
422     def _compute_moments_from_axes(a,b,theta):
423         c = numpy.cos(2.0*theta*ANGLE_UNITS)
424         cp = 1.0+c
425         cm = 1.0-c
426         a_2 = a**2
427         b_2 = b**2
428         return (0.5*(a_2*cp + b_2*cm),0.5*(a_2*cm + b_2*cp),0.5*(a_2-b_2)*numpy.sin(2.0*theta*ANGLE_UNITS))
429
430     @staticmethod
431     def _compute_moments_from_ellipticity(e1,e2,r):
432         e_mag = e1**2 + e2**2
433         half_xx_p_yy = r**2 / (1.0-e_mag)**0.5
434         half_xx_m_yy = e1 * half_xx_p_yy
435         return (half_xx_p_yy+half_xx_m_yy, half_xx_p_yy-half_xx_m_yy, e2*half_xx_p_yy)
436
437     @staticmethod
438     def _compute_axes_from_moments(xx,yy,xy):
439         t = numpy.sqrt((xx-yy)**2+4*xy**2)
440         return (numpy.sqrt(0.5*(xx+yy+t)), numpy.sqrt(0.5*(xx+yy-t)),
441                 0.5*numpy.arctan2(2.0*xy,xx-yy)/ANGLE_UNITS)
442
443     @staticmethod
444     def _compute_axes_from_ellipticity(e1,e2,r):
445         e_mag = (e1**2 + e2**2)**0.5
446         sqrt_q = ((1-e_mag)/(1+e_mag))**0.25
447         return (r/sqrt_q, r*sqrt_q, 0.5*numpy.arctan2(e2,e1)/ANGLE_UNITS)
448
449     @staticmethod
450     def _compute_ellipticity_from_moments(xx,yy,xy):
451         t = xx+yy
452         return (xx-yy)/t, 2*xy/t, (xx*yy-xy**2)**0.25
453
454     @staticmethod
455     def _compute_ellipticity_from_axes(a,b,theta):
456         a_2 = a**2
457         b_2 = b**2
458         m = (a_2-b_2)/(a_2+b_2)
459         return m*numpy.cos(2*theta*ANGLE_UNITS), m*numpy.sin(2*theta*ANGLE_UNITS), numpy.sqrt(a*b)
460