improved masking, now one can add or subtract with GUI; implemented also polygon

xray/mask.py

@@ -5,53 +5,91 @@ import logging
 log = logging.getLogger(__name__)
 
 import os
+import collections
 import numpy as np
 import matplotlib.pyplot as plt
+from matplotlib.path import Path
+
+maskComponent = collections.namedtuple('maskComponent',['operation','geometry','vertices'])
+
+def _rectangleToMask(X,Y,vertices):
+    ( (x1,y1), (x2,y2) ) = vertices
+    if x1>x2: x1,x2=x2,x1
+    if y1>y2: y1,y2=y2,y1
+    return (X>x1) & (X<x2) & ( Y>y1) & (Y<y2)
+
+def _circleToMask(X,Y,vertices):
+    c,p = vertices
+    r = np.sqrt( (p[0]-c[0])**2 + (p[1]-c[1])**2 )
+    d = np.sqrt((X-c[0])**2+(Y-c[1])**2)
+    return d<r
+
+def _polygonToMask(X,Y,vertices):
+    points = np.vstack((X.flatten(),Y.flatten())).T
+    path = Path(vertices)
+    grid = path.contains_points(points)
+    return grid.reshape(X.shape)
 
 class MyMask(object):
   def __init__(self,img=None):
     self.comp = []
     self.img  = img
     self.mask = None
+    self._cache = None
+
+  def _define_component(self,operation,geometry,*vertices):
+    #print("define_comp",type(vertices),vertices)
+    if geometry == 'circle' and len(vertices) == 3:
+      xcen,ycen,radius = vertices
+      vertices = ( (xcen,ycen), (xcen+radius,ycen) )
+    if geometry == 'rectangle' and len(vertices) == 4:
+      vertices = ( (vertices[0],vertices[1]),(vertices[2],vertices[3]) )
+    # make sure vertices tuples
+    if isinstance(vertices,list):
+      vertices = [ (v[0],v[1]) for v in vertices ]
+      vertices = tuple(vertices)
+    a = dict( vertices = None )
+    self.comp.append( maskComponent(operation=operation,vertices=vertices,geometry=geometry) )
+    
+  
+  def addCircle(self,*vertices): self._define_component( 'add', 'circle', *vertices )
+  def subtractCircle(self,*vertices): self._define_component( 'subtract', 'circle', *vertices )
+
+  def addRectangle(self,*vertices): self._define_component( 'add','rectangle', *vertices)
+  def subtractRectangle(self,*vertices): self._define_component( 'subtract','rectangle',*vertices)
+  
+ 
+  def addPolygon(self,*vertices): self._define_component( 'add','polygon',*vertices)
+  def subtractPolygon(self,*vertices): self._define_component( 'subtract','polygon',*vertices)
 
-  def addCircle(self,xcen,ycen,radius):
-    self.comp.append( ["add","circle", [xcen,ycen,radius] ] )
-  def subtractCircle(self,xcen,ycen,radius):
-    self.comp.append( ["subtract","circle", [xcen,ycen,radius] ] )
-  def addRectangle(self,x1,y1,x2,y2):
-    if x1>x2: x1,x2=x2,x1
-    if y1>y2: y1,y2=y2,y1
-    self.comp.append( ["add","rectangle", [x1,y1,x2,y2] ])
-  def subtractRectangle(self,x1,y1,x2,y2):
-    if x1>x2: x1,x2=x2,x1
-    if y1>y2: y1,y2=y2,y1
-    self.comp.append( ["subtract","rectangle", [x1,y1,x2,y2] ])
 
   def getMask(self,shape=None):
-    if shape is None: shape = self.img.shape
+    if shape is None and self.img is not None: shape = self.img.shape
-    m = []
+    if shape is None and self.img is     None: shape = self._cache['shape']
-    X,Y = np.meshgrid ( range(shape[0]),range(shape[1]) )
+    
-    for o in self.comp:
+    if self._cache is None: self._cache = dict( shape = shape )
-      whattodo = o[0]
+    # reset cache if shape does not match
-      kind=o[1]
+    if shape != self._cache['shape']:
-      pars=o[2]
+      self._cache = dict( shape = shape )
-      if kind == "circle":
+    X,Y = np.meshgrid ( range(shape[1]),range(shape[0]) )
-        (xc,yc,r) = pars
+    for component in self.comp:
-        d = np.sqrt((X-xc)**2+(Y-yc)**2)
+      if component not in self._cache:
-        #plt.imshow(d<r)
+        if component.geometry == 'circle':
-        #raw_input()
+          mask = _circleToMask( X,Y,component.vertices )
-        m.append( d<r )
+        elif component.geometry == 'rectangle':
-      if kind == "rectangle":
+          mask = _rectangleToMask( X,Y,component.vertices )
-        (x1,y1,x2,y2) = pars
+        elif component.geometry == 'polygon':
-        temp = (X>x1) & (X<x2) & ( Y>y1) & (Y<y2)
+          mask = _polygonToMask( X,Y,component.vertices )
-        m.append( temp )
-    mask = np.zeros(shape,dtype=np.bool)
-    for i in range(len(m)):
-      whattodo = self.comp[i][0]
-      if (whattodo == "add"):
-        mask[m[i]] = True
         else:
-        mask[m[i]] = False
+          raise ValueError("Mask type %s not recongnized"%component.geometry)
+        self._cache[component] = mask
+    mask = np.zeros(shape,dtype=np.bool)
+    for comp in self.comp:
+      m = self._cache[ comp ]
+      if (comp.operation == "add"):
+        mask[m] = True
+      else:
+        mask[m] = False
     self.mask = mask
     return mask
 
@@ -70,25 +108,19 @@ class MyMask(object):
     i=fabio.edfimage.edfimage(mask.astype(np.uint8)); # edf does not support bool
     i.save(fname)
 
-def test():
-  mask = MyMask()
-  mask.addCircle(400,300,250)
-  mask.subtractCircle(400,300,150)
-  mask.addRectangle(350,250,1500,700)
-  mask.show()
-  return mask
-
 def snap(point,shape,snapRange=20):
   snapped = list(point)
   if snapped[0] < snapRange: snapped[0] = 0
-  if snapped[0] > shape[0]-snapRange: snapped[0] = shape[0]
+  if snapped[0] > shape[1]-snapRange: snapped[0] = shape[1]
   if snapped[1] < snapRange: snapped[1] = 0
-  if snapped[1] > shape[1]-snapRange: snapped[1] = shape[1]
+  if snapped[1] > shape[0]-snapRange: snapped[1] = shape[0]
-  return snapped
+  return tuple(snapped)
 
-def getPoint(shape,snapRange):
+def getPoints(N=1,shape=(100,100),snapRange=0):
-  c = plt.ginput()[0]
+  if N<1: print('Right click cancels last point, middle click ends the polygon')
-  c = snap(c,shape,snapRange=snapRange)
+  c = plt.ginput(N)
+  c = [ snap(point,shape,snapRange=snapRange) for point in c ]
+  if len(c) == 1: c = c[0]
   return c
 
 def makeMaskGui(img,snapRange=60):
@@ -100,20 +132,32 @@ def makeMaskGui(img,snapRange=60):
     plt.clim(np.percentile(img,(2,98)))
     plt.imshow(mask.getMatplotlibMask())
     plt.pause(0.01)
-    ans = input("What's next c/r/done? ")
+    ans = input("What's next p/P/c/C/r/R/done? (capitals = subtract)")
     if ans == "c":
-      print("Adding circle, click on center")
+      print("Adding circle, click on center then another point to define radius")
-      c = getPoint(img.shape,snapRange)
+      vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
-      print("Adding circle, click on another point to define radius")
+      mask.addCircle(*vertices)
-      p = getPoint(img.shape,snapRange)
+    if ans == "C":
-      r = np.sqrt( (p[0]-c[0])**2 + (p[1]-c[1])**2 )
+      print("Subtracting circle, click on center then another point to define radius")
-      mask.addCircle(c[0],c[1],r)
+      vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
+      mask.subtractCircle(*vertices)
     if ans == "r":
-      print("Adding rectangle, click on one corner")
+      print("Adding rectangle, click on one corner and then on the opposite one")
-      c1 = getPoint(img.shape,snapRange)
+      vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
-      print("Adding rectangle, click on opposite corner")
+      mask.addRectangle(*vertices)
-      c2 = getPoint(img.shape,snapRange)
+    if ans == "R":
-      mask.addRectangle(c1[0],c1[1],c2[0],c2[1])
+      print("Subtracting rectangle, click on one corner and then on the opposite one")
+      vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
+      mask.subtractRectangle(*vertices)
+    if ans == 'p':
+      print("Adding polygon")
+      vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange)
+      mask.addPolygon(*vertices)
+    if ans == 'P':
+      print("Subtracting polygon")
+      vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange)
+      mask.subtractPolygon(*vertices)
+
     plt.imshow(mask.getMatplotlibMask())
     plt.pause(0.01)
   fname = input("Enter a valid filename (ext .edf or .npy) if you want to save the mask (empty otherwise)")
@@ -150,6 +194,16 @@ def maskCenterLines(width,shape):
   return mask
 
 
+def test(shape=(1000,2000)):
+  mask = MyMask()
+  mask.addCircle(400,300,250)
+  mask.subtractCircle(400,300,150)
+  mask.addRectangle(350,250,1500,700)
+  plt.imshow( mask.getMask(shape) )
+  return mask
+
+
+
 if __name__ == "__main__":
   test()
   plt.show()