fill_between_posneg.py

#!/usr/bin/env python
"""
From: James Boyle <boyle5@llnl.gov>
Subject: possible candidate for examples directory using fill
To: John Hunter <jdhunter@nitace.bsd.uchicago.edu>
Date: Tue, 8 Mar 2005 15:44:11 -0800

I often compare the output from two sensors and I am interested in the
sign of the differences.

I wrote the enclosed code for find the polygons of the positive and
negative regions of the difference of two curves. It is written for
simple-minded straightforwardness rather than speed or elegance.
It is easy to fill in the two sets of polygons with contrasting colors.
For efficiency one could use collections but my curves are such that
fill is quick enough.

The code uses a simple linear technique to find the crossover point,
this too could be made more sophisticated if one desired.

I have found this code to be very handy for the comparisons I perform
-
maybe someone else would find it useful.

--Jim
"""

#!/usr/bin/env python

from pylab import *

def findZero(i,x,y1,y2):
     im1 = i-1
     m1 = (y1[i] - y1[im1])/(x[i] - x[im1])
     m2 = (y2[i] - y2[im1])/(x[i] - x[im1])
     b1 = y1[im1] - m1*x[im1]
     b2 = y2[im1] - m2*x[im1]
     xZero = (b1 - b2)/(m2 - m1)
     yZero = m1*xZero + b1
     return (xZero, yZero)

def posNegFill(x,y1,y2):
      diff = y2 - y1
      pos = []
      neg = []
      xx1 = [x[0]]
      xx2 = [x[0]]
      yy1 = [y1[0]]
      yy2 = [y2[0]]
      oldSign = (diff[0] < 0 )
      npts = x.shape[0]
      for i in range(1,npts):
          newSign = (diff[i] < 0)
          if newSign != oldSign:
              xz,yz = findZero(i,x,y1,y2)
              xx1.append(xz)
              yy1.append(yz)
              xx2.reverse()
              xx1.extend(xx2)
              yy2.reverse()
              yy1.extend(yy2)
              if oldSign:
                  neg.append( (xx1,yy1) )
              else:
                  pos.append( (xx1,yy1) )
              xx1 = [xz,x[i]]
              xx2 = [xz,x[i]]
              yy1 = [yz,y1[i]]
              yy2 = [yz,y2[i]]
              oldSign = newSign
          else:
              xx1.append( x[i])
              xx2.append( x[i])
              yy1.append(y1[i])
              yy2.append(y2[i])
              if i == npts-1:
                  xx2.reverse()
                  xx1.extend(xx2)
                  yy2.reverse()
                  yy1.extend(yy2)
                  if oldSign :
                      neg.append( (xx1,yy1) )
                  else:
                      pos.append( (xx1,yy1) )
      return pos,neg

x1 = arange(0, 2, 0.01)
y1 = sin(2*pi*x1)
y2 = sin(4*pi*x1)

# find positive and negative polygons of difference
pos,neg =  posNegFill(x1,y1,y2)
# positive y2 > y1 is blue
for x,y  in pos:
     p = fill(x,y,'b')

# negative Y2 < y1 is red
for x,y in neg:
     p = fill(x,y,'r')

show()