diff --git a/Gpyr.py b/Gpyr.py
index 0db48f3..9408e74 100644
--- a/Gpyr.py
+++ b/Gpyr.py
@@ -1,8 +1,8 @@
-from Lpyr import Lpyr
-from namedFilter import namedFilter
-from maxPyrHt import maxPyrHt
+from .Lpyr import Lpyr
+from .namedFilter import namedFilter
+from .maxPyrHt import maxPyrHt
 import numpy
-from corrDn import corrDn
+from .corrDn import corrDn
 
 class Gpyr(Lpyr):
     filt = ''
@@ -13,8 +13,8 @@ class Gpyr(Lpyr):
     def __init__(self, *args):    # (image, height, filter, edges)
         self.pyrType = 'Gaussian'
         if len(args) < 1:
-            print "pyr = Gpyr(image, height, filter, edges)"
-            print "First argument (image) is required"
+            print("pyr = Gpyr(image, height, filter, edges)")
+            print("First argument (image) is required")
             return
         else:
             self.image = args[0]
@@ -22,10 +22,10 @@ def __init__(self, *args):    # (image, height, filter, edges)
         if len(args) > 2:
             filt = args[2]
             if not (filt.shape == 1).any():
-                print "Error: filt should be a 1D filter (i.e., a vector)"
+                print("Error: filt should be a 1D filter (i.e., a vector)")
                 return
         else:
-            print "no filter set, so filter is binom5"
+            print("no filter set, so filter is binom5")
             filt = namedFilter('binom5')
             if self.image.shape[0] == 1:
                 filt = filt.reshape(1,5)
@@ -40,8 +40,8 @@ def __init__(self, *args):    # (image, height, filter, edges)
             else:
                 self.height = args[1]
                 if self.height > maxHeight:
-                    print ( "Error: cannot build pyramid higher than %d levels"
-                            % (maxHeight) )
+                    print(( "Error: cannot build pyramid higher than %d levels"
+                            % (maxHeight) ))
                     return
         else:
             self.height = maxHeight
diff --git a/JBhelpers.py b/JBhelpers.py
index 9261b77..e097382 100644
--- a/JBhelpers.py
+++ b/JBhelpers.py
@@ -1,6 +1,6 @@
 class struct( object ):
     def __init__( self, **kwargs ):
-        for k, v in kwargs.iteritems():
+        for k, v in kwargs.items():
             setattr( self, k, v )
 
 
@@ -65,7 +65,7 @@ def nbimage( data, vmin = None, vmax = None, vsym = False, saveas = None ):
     '''
     from IPython.display import display, Image
     from PIL.Image import fromarray
-    from StringIO import StringIO
+    from io import StringIO
     data = rerange( data, vmin, vmax, vsym )
     data = data.squeeze()
     # try to be smart
@@ -93,14 +93,14 @@ def nbimageLCVbak( data, vmin = None, vmax = None, vsym = False, saveas = None,
     '''
     from IPython.display import display, Image, HTML
     from PIL.Image import fromarray
-    from StringIO import StringIO
+    from io import StringIO
     css_styling()
     data = rerange( data, vmin, vmax, vsym )
     data = data.squeeze()
 
     # try to be smart
     if 3 <= data.shape[ 0 ] <= 4:
-        print 'transposing'
+        print('transposing')
         data = data.transpose( ( 1, 2, 0 ) )
 
     s = StringIO()
@@ -135,7 +135,7 @@ def nbimageLCVbak2( data, vmin = None, vmax = None, vsym = False, saveas = None,
     '''
     from IPython.display import display, Image, HTML
     from PIL.Image import fromarray
-    from StringIO import StringIO
+    from io import StringIO
     import base64
     from PyQt4 import QtGui
     from PyQt4 import QtCore
@@ -145,7 +145,7 @@ def nbimageLCVbak2( data, vmin = None, vmax = None, vsym = False, saveas = None,
     data = data.squeeze()
     # try to be smart
     if 3 <= data.shape[ 0 ] <= 4:
-        print 'transposing'
+        print('transposing')
         data = data.transpose( ( 1, 2, 0 ) )
     s = StringIO()
     fromarray( data ).save( s, 'png' )
@@ -157,13 +157,13 @@ def nbimageLCVbak2( data, vmin = None, vmax = None, vsym = False, saveas = None,
 
     matrix = numpy.require(data, numpy.uint8, 'C')
     (w, h) = matrix.shape
-    print matrix
+    print(matrix)
     qim = QtGui.QImage(matrix.data, w, h, QtGui.QImage.Format_Indexed8)
     qim.ndarray = matrix    # do we need this?
     
     # make colormap
     incr = (256/nshades)+1
-    colors = range(0,255,(256/nshades)+1)
+    colors = list(range(0,255,(256/nshades)+1))
     colors[-1] = 255
     colctr = -1
     for i in range(256):
@@ -216,7 +216,7 @@ def nbimageLCV( dlist, vmin = None, vmax = None, vsym = False, saveas = None,
     '''
     from IPython.display import display, Image, HTML
     from PIL.Image import fromarray
-    from StringIO import StringIO
+    from io import StringIO
     import base64
     from PyQt4 import QtGui
     from PyQt4 import QtCore
@@ -255,7 +255,7 @@ def nbimageLCV( dlist, vmin = None, vmax = None, vsym = False, saveas = None,
     
         # make colormap
         incr = (256/nshades)+1
-        colors = range(0,255,(256/nshades)+1)
+        colors = list(range(0,255,(256/nshades)+1))
         colors[-1] = 255
         colctr = -1
         for i in range(256):
@@ -325,7 +325,7 @@ def showIm( dlist, v = None, zoom = 1, title = "", nshades = 256, ncols = 1):
     '''
     from IPython.display import display, Image, HTML
     from PIL.Image import fromarray
-    from StringIO import StringIO
+    from io import StringIO
     import base64
     from PyQt4 import QtGui
     from PyQt4 import QtCore
@@ -353,9 +353,9 @@ def showIm( dlist, v = None, zoom = 1, title = "", nshades = 256, ncols = 1):
             vmin = p1-(p2-p1)/8.0
             vmax = p2+(p2-p1)/8.0
         else:
-            print "Error: range of %s is not recognized." % v
-            print "       please use a two element tuple or "
-            print "       'auto', 'auto2' or 'auto3'"
+            print("Error: range of %s is not recognized." % v)
+            print("       please use a two element tuple or ")
+            print("       'auto', 'auto2' or 'auto3'")
             return
 
         data = rerange( data, vmin, vmax, vsym )
@@ -377,7 +377,7 @@ def showIm( dlist, v = None, zoom = 1, title = "", nshades = 256, ncols = 1):
     
         # make colormap
         incr = (256/nshades)+1
-        colors = range(0,255,(256/nshades)+1)
+        colors = list(range(0,255,(256/nshades)+1))
         colors[-1] = 255
         colctr = -1
         for i in range(256):
diff --git a/Lpyr.py b/Lpyr.py
index d8992b7..1773d7a 100644
--- a/Lpyr.py
+++ b/Lpyr.py
@@ -1,10 +1,10 @@
-from pyramid import pyramid
-from corrDn import corrDn
-from namedFilter import namedFilter
-from maxPyrHt import maxPyrHt
-from upConv import upConv
-from showIm import showIm
-import JBhelpers
+from .pyramid import pyramid
+from .corrDn import corrDn
+from .namedFilter import namedFilter
+from .maxPyrHt import maxPyrHt
+from .upConv import upConv
+from .showIm import showIm
+from . import JBhelpers
 import numpy
 import math
 import matplotlib
@@ -20,17 +20,17 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
         if len(args) > 0:
             self.image = args[0]
         else:
-            print "pyr = Lpyr(image, height, filter1, filter2, edges)"
-            print "First argument (image) is required"
+            print("pyr = Lpyr(image, height, filter1, filter2, edges)")
+            print("First argument (image) is required")
             return
 
         if len(args) > 2:
             filt1 = args[2]
-            if isinstance(filt1, basestring):
+            if isinstance(filt1, str):
                 filt1 = namedFilter(filt1)
             elif len(filt1.shape) != 1 and ( filt1.shape[0] != 1 and
                                              filt1.shape[1] != 1 ):
-                print "Error: filter1 should be a 1D filter (i.e., a vector)"
+                print("Error: filter1 should be a 1D filter (i.e., a vector)")
                 return
         else:
             filt1 = namedFilter('binom5')
@@ -41,11 +41,11 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
 
         if len(args) > 3:
             filt2 = args[3]
-            if isinstance(filt2, basestring):
+            if isinstance(filt2, str):
                 filt2 = namedFilter(filt2)
             elif len(filt2.shape) != 1 and ( filt2.shape[0] != 1 and
                                             filt2.shape[1] != 1 ):
-                print "Error: filter2 should be a 1D filter (i.e., a vector)"
+                print("Error: filter2 should be a 1D filter (i.e., a vector)")
                 return
         else:
             filt2 = filt1
@@ -58,8 +58,8 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
             else:
                 self.height = args[1]
                 if self.height > maxHeight:
-                    print ( "Error: cannot build pyramid higher than %d levels"
-                            % (maxHeight) )
+                    print(( "Error: cannot build pyramid higher than %d levels"
+                            % (maxHeight) ))
                     return
         else:
             self.height = maxHeight
@@ -144,8 +144,8 @@ def catBands(self, *args):
     # set a pyramid value
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, element(tuple), value)'
+            print('Error: three input parameters required:')
+            print('  set(band, element(tuple), value)')
         self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
 
     def reconPyr(self, *args):
@@ -169,15 +169,15 @@ def reconPyr(self, *args):
 
         maxLev = self.height
 
-        if isinstance(levs, (str,basestring)) and levs == 'all':
-            levs = range(0,maxLev)
+        if isinstance(levs, str) and levs == 'all':
+            levs = list(range(0,maxLev))
         else:
             if (levs > maxLev-1).any():
-                print ( "Error: level numbers must be in the range [0, %d]." % 
-                        (maxLev-1) )
+                print(( "Error: level numbers must be in the range [0, %d]." % 
+                        (maxLev-1) ))
                 return
 
-        if isinstance(filt2, basestring):
+        if isinstance(filt2, str):
             filt2 = namedFilter(filt2)
         else:
             if len(filt2.shape) == 1:
@@ -297,11 +297,11 @@ def showPyr(self, pRange = None, gap = 1, scale = None, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.std(band)
             pRange[nind,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(pRange, basestring):
-            print "Error: band range argument: %s" % (pRange)
+        elif isinstance(pRange, str):
+            print("Error: band range argument: %s" % (pRange))
             return
         elif pRange.shape[0] == 1 and pRange.shape[1] == 2:
-            scales = numpy.power( numpy.array( range(0,nind) ), scale)
+            scales = numpy.power( numpy.array( list(range(0,nind)) ), scale)
             pRange = numpy.outer( scales, pRange )
             band = self.pyrLow()
             pRange[nind,:] = ( pRange[nind,:] + numpy.mean(band) - 
@@ -349,7 +349,7 @@ def showPyr(self, pRange = None, gap = 1, scale = None, disp = 'qt'):
                 llpos[bnum,:] = ctr - numpy.floor(numpy.array(sz))/2.0 
             # make position list positive, and allocate appropriate image
             llpos = llpos - numpy.ones((nind,1))*numpy.min(llpos)
-            pind = range(self.height)
+            pind = list(range(self.height))
             for i in pind:
                 pind[i] = self.band(i).shape
             urpos = llpos + pind
diff --git a/SCFpyr.py b/SCFpyr.py
index d368cb2..2afee52 100644
--- a/SCFpyr.py
+++ b/SCFpyr.py
@@ -1,211 +1,298 @@
-from SFpyr import SFpyr
-import numpy
-from steer2HarmMtx import steer2HarmMtx
-from rcosFn import rcosFn
-from pointOp import pointOp
+from .SFpyr import SFpyr
+import numpy as np
+from .rcosFn import rcosFn
+from .pointOp import pointOp
 import scipy
-from mkAngle import mkAngle
-import cmath
+from .mkAngle import mkAngle
+
 
 class SCFpyr(SFpyr):
     filt = ''
     edges = ''
     
     #constructor
-    def __init__(self, *args):    # (image, height, order, twidth)
+    def __init__(self, image, height, order, twidth, scale, n_scales, xp=np):    # (image, height, order, twidth, scale, n_scales)
         self.pyrType = 'steerableFrequency'
-
-        if len(args) > 0:
-            self.image = args[0]
-        else:
-            print "First argument (image) is required."
-            return
-
-        #------------------------------------------------
-        # defaults:
-
-        max_ht = numpy.floor( numpy.log2( min(self.image.shape) ) ) - 2
-        if len(args) > 1:
-            if(args[1] > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
-            ht = args[1]
-        else:
-            ht = max_ht
-        ht = int(ht)
-            
-        if len(args) > 2:
-            if args[2] > 15 or args[2] < 0:
-                print "Warning: order must be an integer in the range [0,15]. Truncating."
-                order = min( max(args[2],0), 15 )
-            else:
-                order = args[2]
-        else:
-            order = 3
-
-        nbands = order+1
-
-        if len(args) > 3:
-            if args[3] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
-                twidth = 1
-            else:
-                twidth = args[3]
-        else:
-            twidth = 1
-
-        #------------------------------------------------------
-        # steering stuff:
-
-        if nbands % 2 == 0:
-            harmonics = numpy.array(range(nbands/2)) * 2 + 1
-        else:
-            harmonics = numpy.array(range((nbands-1)/2)) * 2
-
-        steermtx = steer2HarmMtx(harmonics,
-                                 numpy.pi*numpy.array(range(nbands))/nbands,
-                                 'even')
-        #------------------------------------------------------
-        
-        dims = numpy.array(self.image.shape)
-        ctr = numpy.ceil((numpy.array(dims)+0.5)/2)
-        
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
-                                     (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
-                                     (dims[0]/2))
-        angle = numpy.arctan2(yramp, xramp)
-        log_rad = numpy.sqrt(xramp**2 + yramp**2)
+        self.image = image
+        self.ht = height
+        self.order = order
+        self.twidth = twidth
+        self.scale = scale
+        self.n_scales = n_scales
+        self.nbands = self.order+1
+
+        self.xp = xp
+        dims = np.array(self.image.shape)
+        ctr = np.ceil((dims+0.5)/2).astype('int')
+
+        (xramp, yramp) = self.xp.meshgrid((self.xp.arange(1, dims[1]+1)-ctr[1]) / (dims[1]/2.),
+                                          (self.xp.arange(1, dims[0]+1)-ctr[0]) / (dims[0]/2.))
+        angle = self.xp.arctan2(yramp, xramp)
+        log_rad = self.xp.sqrt(xramp**2 + yramp**2)
         log_rad[ctr[0]-1, ctr[1]-1] = log_rad[ctr[0]-1, ctr[1]-2]
-        log_rad = numpy.log2(log_rad);
+        log_rad = self.xp.log2(log_rad)
 
         ## Radial transition function (a raised cosine in log-frequency):
-        (Xrcos, Yrcos) = rcosFn(twidth, (-twidth/2.0), numpy.array([0,1]))
-        Yrcos = numpy.sqrt(Yrcos)
-
-        YIrcos = numpy.sqrt(1.0 - Yrcos**2)
-        lo0mask = pointOp(log_rad, YIrcos, Xrcos[0], Xrcos[1]-Xrcos[0], 0)
+        (Xrcos, Yrcos) = rcosFn(self.twidth, (-self.twidth/2.0), np.array([0,1]))
+        Yrcos = self.xp.sqrt(self.xp.array(Yrcos))
 
-        imdft = numpy.fft.fftshift(numpy.fft.fft2(self.image))
+        YIrcos = self.xp.sqrt(1.0 - Yrcos**2)
 
+        imdft = self.xp.fft.fftshift(self.xp.fft.fft2(self.image))
+        if self.xp.__name__ == 'numpy':
+            lo0mask = pointOp(log_rad, YIrcos, Xrcos[0], Xrcos[1]-Xrcos[0], 0)
+            hi0mask = pointOp(log_rad, Yrcos, Xrcos[0], Xrcos[1] - Xrcos[0], 0)
+        else:
+            lo0mask = self.xp.array(pointOp(self.xp.asnumpy(log_rad), self.xp.asnumpy(YIrcos), self.xp.asnumpy(Xrcos)[0],
+                                            self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
+            hi0mask = self.xp.array(pointOp(self.xp.asnumpy(log_rad), self.xp.asnumpy(Yrcos), self.xp.asnumpy(Xrcos)[0],
+                                            self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
         self.pyr = []
         self.pyrSize = []
 
-        hi0mask = pointOp(log_rad, Yrcos, Xrcos[0], Xrcos[1]-Xrcos[0], 0)
-
         hi0dft = imdft * hi0mask.reshape(imdft.shape[0], imdft.shape[1])
-        hi0 = numpy.fft.ifft2(numpy.fft.ifftshift(hi0dft))
+        hi0 = self.xp.fft.ifft2(self.xp.fft.ifftshift(hi0dft))
 
-        self.pyr.append(numpy.real(hi0.copy()))
+        self.pyr.append(self.xp.real(hi0.copy()))
         self.pyrSize.append(hi0.shape)
 
         lo0mask = lo0mask.reshape(imdft.shape[0], imdft.shape[1])
         lodft = imdft * lo0mask
 
-        for i in range(ht):
-            bands = numpy.zeros((lodft.shape[0]*lodft.shape[1], nbands))
-            bind = numpy.zeros((nbands, 2))
-        
-            Xrcos -= numpy.log2(2)
+        # self.pind = numpy.zeros((nbands, 2))
+        self.bands = []
+        for i in range(self.ht - 1, -1, -1):
+            # Xrcos -= numpy.log2(2)
+            Xrcos -= np.log2(1. / self.scale)
 
             lutsize = 1024
-            Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
-
-            order = nbands -1
-            const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
-
-            alfa = ( (numpy.pi+Xcosn) % (2.0*numpy.pi) ) - numpy.pi
-            Ycosn = ( 2.0*numpy.sqrt(const) * (numpy.cos(Xcosn)**order) * 
-                      (numpy.abs(alfa)<numpy.pi/2.0).astype(int) )
-            log_rad_tmp = numpy.reshape(log_rad, (1,log_rad.shape[0]*
-                                                  log_rad.shape[1]))
-            himask = pointOp(log_rad_tmp, Yrcos, Xrcos[0], Xrcos[1]-Xrcos[0], 0)
-            
+            Xcosn = np.pi * self.xp.array(self.xp.arange(-(2*lutsize+1), (lutsize+2))) / lutsize
+
+            order = self.nbands - 1
+            const = (2**(2*order))*(scipy.special.factorial(order, exact=True)**2) / \
+                    float(self.nbands*scipy.special.factorial(2*order, exact=True))
+
+            alfa = ( (np.pi+Xcosn) % (2.0*np.pi) ) - np.pi
+            Ycosn = 2.0 * self.xp.sqrt(const) * (self.xp.cos(Xcosn)**order) * (self.xp.abs(alfa) < np.pi/2.0).astype(int)
+            log_rad_tmp = self.xp.reshape(log_rad, (1,log_rad.shape[0]* log_rad.shape[1]))
+
+            if self.xp.__name__ == 'numpy':
+                himask = pointOp(log_rad_tmp, Yrcos, Xrcos[0], Xrcos[1]-Xrcos[0], 0)
+            else:
+                himask = self.xp.array(pointOp(self.xp.asnumpy(log_rad_tmp), self.xp.asnumpy(Yrcos), self.xp.asnumpy(Xrcos)[0],
+                                               self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
+
             himask = himask.reshape(lodft.shape[0], lodft.shape[1])
-            for b in range(nbands):
-                angle_tmp = numpy.reshape(angle, 
-                                          (1,angle.shape[0]*angle.shape[1]))
-                anglemask = pointOp(angle_tmp, Ycosn,
-                                    Xcosn[0]+numpy.pi*b/nbands, 
-                                    Xcosn[1]-Xcosn[0], 0)
+            for b in range(self.nbands):
+                angle_tmp = self.xp.reshape(angle, (1, angle.shape[0]*angle.shape[1]))
+
+                if self.xp.__name__ == 'numpy':
+                    anglemask = pointOp(angle_tmp, Ycosn, Xcosn[0] + np.pi*b/self.nbands, Xcosn[1] - Xcosn[0], 0)
+                else:
+                    anglemask = self.xp.array(pointOp(self.xp.asnumpy(angle_tmp), self.xp.asnumpy(Ycosn),
+                                                      self.xp.asnumpy(Xcosn)[0] + np.pi * b / self.nbands,
+                                                      self.xp.asnumpy(Xcosn)[1] - self.xp.asnumpy(Xcosn)[0], 0))
                 anglemask = anglemask.reshape(lodft.shape[0], lodft.shape[1])
-                banddft = (cmath.sqrt(-1)**order) * lodft * anglemask * himask
-                band = numpy.negative(numpy.fft.ifft2(numpy.fft.ifftshift(banddft)))
+                banddft = (self.xp.sqrt(-1 + 0j)**order) * lodft * anglemask * himask
+                band = self.xp.negative(self.xp.fft.ifft2(self.xp.fft.ifftshift(banddft)))
+                self.bands.append(band[:].reshape(-1))
                 self.pyr.append(band.copy())
                 self.pyrSize.append(band.shape)
-
-            dims = numpy.array(lodft.shape)
-            ctr = numpy.ceil((dims+0.5)/2)
-            lodims = numpy.ceil((dims-0.5)/2)
-            loctr = numpy.ceil((lodims+0.5)/2)
+            dims = np.array(lodft.shape)
+            ctr = np.ceil((dims+0.5)/2).astype('int')
+            lodims = np.round(dims * self.scale ** (self.n_scales - self.ht + 1))
+            loctr = np.ceil((lodims+0.5)/2).astype('int')
             lostart = ctr - loctr
-            loend = lostart + lodims
+            loend = (lostart + lodims).astype('int')
 
             log_rad = log_rad[lostart[0]:loend[0], lostart[1]:loend[1]]
             angle = angle[lostart[0]:loend[0], lostart[1]:loend[1]]
             lodft = lodft[lostart[0]:loend[0], lostart[1]:loend[1]]
-            YIrcos = numpy.abs(numpy.sqrt(1.0 - Yrcos**2))
-            log_rad_tmp = numpy.reshape(log_rad, 
-                                        (1,log_rad.shape[0]*log_rad.shape[1]))
-            lomask = pointOp(log_rad_tmp, YIrcos, Xrcos[0], Xrcos[1]-Xrcos[0],
-                             0)
+            YIrcos = self.xp.abs(self.xp.sqrt(1.0 - Yrcos**2))
+            log_rad_tmp = self.xp.reshape(log_rad, (1,log_rad.shape[0]*log_rad.shape[1]))
+
+            if self.xp.__name__ == 'numpy':
+                lomask = pointOp(log_rad_tmp, YIrcos, Xrcos[0], Xrcos[1]-Xrcos[0],0)
+            else:
+                lomask = self.xp.array(pointOp(self.xp.asnumpy(log_rad_tmp), self.xp.asnumpy(YIrcos), self.xp.asnumpy(Xrcos)[0],
+                                               self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
             lodft = lodft * lomask.reshape(lodft.shape[0], lodft.shape[1])
 
-        lodft = numpy.fft.ifft2(numpy.fft.ifftshift(lodft))
-        self.pyr.append(numpy.real(numpy.array(lodft).copy()))
+        self.bands = self.xp.concatenate(self.bands, 0)
+        lodft = self.xp.fft.ifft2(self.xp.fft.ifftshift(lodft))
+        self.pyr.append(self.xp.real(self.xp.array(lodft).copy()))
         self.pyrSize.append(lodft.shape)
+        self.pyrSize = np.array(self.pyrSize)
 
-    # methods
-    def reconPyr(self, *args):
-        if len(args) > 0:
-            levs = args[0]
-        else:
-            levs = 'all'
+    def band(self, bandNum):
+        return self.xp.array(self.pyr[bandNum])
 
-        if len(args) > 1:
-            bands = args[1]
-        else:
-            bands = 'all'
+    def spyrHt(self):
+        return self.ht
 
-        if len(args) > 2:
-            if args[2] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
-                twidth = 1
-            else:
-                twidth = args[2]
-        else:
-            twidth = 1
+    def numBands(self):
+        return self.order + 1
+
+    def pyrLow(self):
+        return self.xp.array(self.band(len(self.pyrSize)-1))
 
-        #-----------------------------------------------------------------
+    def pyrHigh(self):
+        return self.xp.array(self.band(0))
 
+    # methods
+    def reconPyr(self, levs='all', bands='all'):
         pind = self.pyrSize
-        Nsc = int(numpy.log2(pind[0][0] / pind[-1][0]))
-        Nor = (len(pind)-2) / Nsc
+        Nsc = int(np.round(np.log2(pind[0, 0] / pind[-1, 0].astype('float')) / np.log2(1 / self.scale)))
+        Nor = (len(pind)-2) // Nsc
 
         pyrIdx = 1
         for nsc in range(Nsc):
-            firstBnum = nsc * Nor+2
+            firstBnum = nsc * Nor + 1
             dims = pind[firstBnum][:]
-            ctr = (numpy.ceil((dims[0]+0.5)/2.0), numpy.ceil((dims[1]+0.5)/2.0)) #-1?
-            ang = mkAngle(dims, 0, ctr)
-            ang[ctr[0]-1, ctr[1]-1] = -numpy.pi/2.0
+            ctr = np.ceil((dims + 0.5) / 2).astype('int') #-1?
+            ang = self.xp.array(mkAngle(dims, 0, ctr))
+            ang[ctr[0]-1, ctr[1]-1] = -np.pi / 2.0
             for nor in range(Nor):
                 nband = nsc * Nor + nor + 1
                 ch = self.pyr[nband]
-                ang0 = numpy.pi * nor / Nor
-                xang = ((ang-ang0+numpy.pi) % (2.0*numpy.pi)) - numpy.pi
-                amask = 2 * (numpy.abs(xang) < (numpy.pi/2.0)).astype(int) + (numpy.abs(xang) == (numpy.pi/2.0)).astype(int)
+                ang0 = np.pi * nor / Nor
+                xang = ((ang - ang0 + np.pi) % (2.0 * np.pi)) - np.pi
+                amask = 2 * (self.xp.abs(xang) < (np.pi/2.0)).astype(int) + (self.xp.abs(xang) == (np.pi/2.0)).astype(int)
                 amask[ctr[0]-1, ctr[1]-1] = 1
                 amask[:,0] = 1
                 amask[0,:] = 1
-                amask = numpy.fft.fftshift(amask)
-                ch = numpy.fft.ifft2(amask * numpy.fft.fft2(ch))  # 'Analytic' version
+                amask = self.xp.fft.fftshift(amask)
+                ch = self.xp.fft.ifft2(amask * self.xp.fft.fft2(ch))  # 'Analytic' version
                 # f = 1.000008  # With this factor the reconstruction SNR
                                 # goes up around 6 dB!
                 f = 1
-                ch = f*0.5*numpy.real(ch)   # real part
+                ch = f*0.5*self.xp.real(ch)   # real part
                 self.pyr[pyrIdx] = ch
                 pyrIdx += 1
 
-        res = self.reconSFpyr(levs, bands, twidth);
+        res = self._reconSFpyr(levs, bands)
 
         return res
+
+    def _reconSFpyr(self, levs='all', bands='all'):
+        nbands = self.numBands()
+
+        maxLev = 1 + self.spyrHt()
+        if isinstance(levs, str) and levs == 'all':
+            levs = np.array(list(range(maxLev + 1)))
+        elif isinstance(levs, str):
+            raise ValueError("Error: %s not valid for levs parameter. "
+                             "levs must be either a 1D numpy array or the string 'all'." % (levs))
+        else:
+            levs = np.array(levs)
+
+        if isinstance(bands, str) and bands == 'all':
+            bands = np.array(list(range(nbands)))
+        elif isinstance(bands, str):
+            raise ValueError("Error: %s not valid for bands parameter. "
+                             "bands must be either a 1D numpy array or the string 'all'." % (bands))
+        else:
+            bands = np.array(bands)
+
+        # -------------------------------------------------------------------
+        # matlab code starts here
+        pind = self.pyrSize
+        dims = np.array(self.pyrSize[0])
+        ctr = np.ceil((dims + 0.5) / 2.0).astype('int')
+
+        (xramp, yramp) = self.xp.meshgrid((self.xp.arange(1, dims[1] + 1) - ctr[1]) / (dims[1] / 2.),
+                                        (self.xp.arange(1, dims[0] + 1) - ctr[0]) / (dims[0] / 2.))
+        angle = self.xp.arctan2(yramp, xramp)
+        log_rad = self.xp.sqrt(xramp ** 2 + yramp ** 2)
+        log_rad[ctr[0] - 1, ctr[1] - 1] = log_rad[ctr[0] - 1, ctr[1] - 2]
+        log_rad = self.xp.log2(log_rad)
+
+        ## Radial transition function (a raised cosine in log-frequency):
+        (Xrcos, Yrcos) = rcosFn(self.twidth, (-self.twidth / 2.0), np.array([0, 1]))
+        Xrcos = self.xp.array(Xrcos)
+        Yrcos = self.xp.sqrt(self.xp.array(Yrcos))
+        YIrcos = self.xp.sqrt(self.xp.abs(1.0 - Yrcos ** 2))
+
+        resdft = self._reconvSFpyr(self.pyr[1:], pind[1:], log_rad, Xrcos, Yrcos, angle, nbands, levs, bands, self.scale)
+
+        # apply lo0mask
+        if self.xp.__name__ == 'numpy':
+            lo0mask = pointOp(log_rad, YIrcos, Xrcos[0], Xrcos[1] - Xrcos[0], 0)
+        else:
+            lo0mask = self.xp.array(pointOp(self.xp.asnumpy(log_rad), self.xp.asnumpy(YIrcos), self.xp.asnumpy(Xrcos)[0],
+                                            self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
+        resdft = resdft * lo0mask
+
+        # residual highpass subband
+        if self.xp.__name__ == 'numpy':
+            hi0mask = pointOp(log_rad, Yrcos, Xrcos[0], Xrcos[1] - Xrcos[0], 0)
+        else:
+            hi0mask = self.xp.array(pointOp(self.xp.asnumpy(log_rad), self.xp.asnumpy(Yrcos), self.xp.asnumpy(Xrcos)[0],
+                                            self.xp.asnumpy(Xrcos)[1] - self.xp.asnumpy(Xrcos)[0], 0))
+
+        hi0mask = hi0mask.reshape(resdft.shape[0], resdft.shape[1])
+        if 0 in levs:
+            resdft += self.xp.fft.fftshift(self.xp.fft.fft2(self.pyr[0])) * hi0mask
+
+        outresdft = self.xp.real(self.xp.fft.ifft2(self.xp.fft.ifftshift(resdft)))
+        return outresdft
+
+    def _reconvSFpyr(self, pyr, pind, log_rad, Xrcos, Yrcos, angle, nbands, levs, bands, scale):
+        lo_ind = nbands + 1
+        dims = pind[0]
+        ctr = np.ceil((dims + 0.5) / 2).astype('int')
+        XXrcos = self.xp.copy(Xrcos - self.xp.log2(1. / scale))
+
+        if (levs > 1).any():
+            lodims = pind[bands[-1] + 1]
+            loctr = np.ceil((lodims + 0.5) / 2).astype('int')
+            lostart = ctr - loctr + 1
+            loend = lostart + lodims - 1
+            nlog_rad = self.xp.copy(log_rad[lostart[0] - 1:loend[0], lostart[1] - 1:loend[1]])
+            nangle = self.xp.copy(angle[lostart[0] - 1: loend[0], lostart[1] - 1: loend[1]])
+
+            if pind.shape[0] > lo_ind:
+                nresdft = self._reconvSFpyr(pyr[lo_ind - 1:len(pyr)], pind[lo_ind - 1:], nlog_rad, XXrcos, Yrcos,
+                                            nangle, nbands, levs - 1, bands, scale)
+            else:
+                nresdft = self.xp.fft.fftshift(self.xp.fft.fft2(pyr[lo_ind - 1]))
+
+            YIrcos = self.xp.sqrt(self.xp.abs(1. - Yrcos ** 2))
+
+            if self.xp.__name__ == 'numpy':
+                lomask = pointOp(nlog_rad, YIrcos, XXrcos[0], XXrcos[1] - XXrcos[0], 0)
+            else:
+                lomask = self.xp.array(pointOp(self.xp.asnumpy(nlog_rad), self.xp.asnumpy(YIrcos), self.xp.asnumpy(XXrcos)[0],
+                                               self.xp.asnumpy(XXrcos)[1] - self.xp.asnumpy(XXrcos)[0], 0))
+
+            resdft = self.xp.zeros(dims).astype('complex128')
+            resdft[lostart[0] - 1:loend[0], lostart[1] - 1:loend[1]] = nresdft * lomask
+        else:
+            resdft = self.xp.zeros(dims)
+
+        if (levs == 1).any():
+            lutsize = 1024
+            Xcosn = np.pi * self.xp.arange(-(2*lutsize + 1), lutsize + 1) / lutsize
+            order = nbands - 1
+            const = 2 ** (2*order) * (scipy.special.factorial(order) ** 2) / (nbands * scipy.special.factorial(2.*order))
+            Ycosn = self.xp.sqrt(const) * (self.xp.cos(Xcosn)) ** order
+
+            if self.xp.__name__ == 'numpy':
+                himask = pointOp(log_rad, Yrcos, XXrcos[0], XXrcos[1] - XXrcos[0], 0)
+            else:
+                himask = self.xp.array(pointOp(self.xp.asnumpy(log_rad), self.xp.asnumpy(Yrcos), self.xp.asnumpy(XXrcos)[0],
+                                               self.xp.asnumpy(XXrcos)[1] - self.xp.asnumpy(XXrcos)[0], 0))
+
+            ind = 0
+            for b in range(nbands):
+                if (bands == b).any():
+                    if self.xp.__name__ == 'numpy':
+                        anglemask = pointOp(angle, Ycosn, Xcosn[0] + np.pi*b/nbands, Xcosn[1] - Xcosn[0], 0)
+                    else:
+                        anglemask = self.xp.asarray(pointOp(self.xp.asnumpy(angle), self.xp.asnumpy(Ycosn), self.xp.asnumpy(Xcosn)[0] +
+                                                       np.pi * b / nbands, self.xp.asnumpy(Xcosn)[1] - self.xp.asnumpy(Xcosn)[0], 0))
+                    band = pyr[ind]
+                    banddft = self.xp.fft.fftshift(self.xp.fft.fft2(band))
+                    resdft += (self.xp.power(-1 + 0j, 0.5)) ** (nbands - 1) * banddft * anglemask * himask
+                ind += 1
+        return resdft
diff --git a/SFpyr.py b/SFpyr.py
index d649a47..ed42697 100644
--- a/SFpyr.py
+++ b/SFpyr.py
@@ -1,8 +1,8 @@
-from Spyr import Spyr
+from .Spyr import Spyr
 import numpy
-from steer2HarmMtx import steer2HarmMtx
-from rcosFn import rcosFn
-from pointOp import pointOp
+from .steer2HarmMtx import steer2HarmMtx
+from .rcosFn import rcosFn
+from .pointOp import pointOp
 import scipy
 
 class SFpyr(Spyr):
@@ -16,7 +16,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         if len(args) > 0:
             self.image = args[0]
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -25,7 +25,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         max_ht = numpy.floor( numpy.log2( min(self.image.shape) ) ) - 2
         if len(args) > 1:
             if(args[1] > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (max_ht))
             ht = args[1]
         else:
             ht = max_ht
@@ -33,7 +33,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             
         if len(args) > 2:
             if args[2] > 15 or args[2] < 0:
-                print "Warning: order must be an integer in the range [0,15]. Truncating."
+                print("Warning: order must be an integer in the range [0,15]. Truncating.")
                 order = min( max(args[2],0), 15 )
             else:
                 order = args[2]
@@ -44,7 +44,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
 
         if len(args) > 3:
             if args[3] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
+                print("Warning: twidth must be positive. Setting to 1.")
                 twidth = 1
             else:
                 twidth = args[3]
@@ -55,21 +55,21 @@ def __init__(self, *args):    # (image, height, order, twidth)
         # steering stuff:
 
         if nbands % 2 == 0:
-            harmonics = numpy.array(range(nbands/2)) * 2 + 1
+            harmonics = numpy.array(list(range(nbands/2))) * 2 + 1
         else:
-            harmonics = numpy.array(range((nbands-1)/2)) * 2
+            harmonics = numpy.array(list(range((nbands-1)/2))) * 2
 
         steermtx = steer2HarmMtx(harmonics, numpy.pi *
-                                 numpy.array(range(nbands))/nbands, 'even')
+                                 numpy.array(list(range(nbands)))/nbands, 'even')
 
         #------------------------------------------------------
         
         dims = numpy.array(self.image.shape)
         ctr = numpy.ceil((numpy.array(dims)+0.5)/2)
         
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
+        (xramp, yramp) = numpy.meshgrid((numpy.array(list(range(1,dims[1]+1)))-ctr[1])/
                                      (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
+                                     (numpy.array(list(range(1,dims[0]+1)))-ctr[0])/
                                      (dims[0]/2))
         angle = numpy.arctan2(yramp, xramp)
         log_rad = numpy.sqrt(xramp**2 + yramp**2)
@@ -106,7 +106,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             Xrcos -= numpy.log2(2)
 
             lutsize = 1024
-            Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
+            Xcosn = numpy.pi * numpy.array(list(range(-(2*lutsize+1), (lutsize+2)))) / lutsize
 
             order = nbands -1
             const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
@@ -160,8 +160,8 @@ def numBands(self):      # FIX: why isn't this inherited
             return 0
         else:
             b = 2
-            while ( b <= len(self.pyrSize) and 
-                    self.pyrSize[b] == self.pyrSize[1] ):
+            while ( b < len(self.pyrSize) and
+                        (self.pyrSize[b] == self.pyrSize[1]).all() ):
                 b += 1
             return b-1
 
@@ -172,47 +172,26 @@ def spyrHt(self):
             spHt = 0
         return spHt
 
-    def reconSFpyr(self, *args):
-        
-        if len(args) > 0:
-            levs = args[0]
-        else:
-            levs = 'all'
-
-        if len(args) > 1:
-            bands = args[1]
-        else:
-            bands = 'all'
-
-        if len(args) > 2:
-            if args[2] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
-                twidth = 1
-            else:
-                twidth = args[2]
-        else:
-            twidth = 1
-
+    def reconSFpyr(self, levs='all', bands='all', twidth=1):
+        assert twidth > 0, 'twidth must be a positive number.'
         #-----------------------------------------------------------------
 
         nbands = self.numBands()
         
         maxLev = 1 + self.spyrHt()
-        if isinstance(levs, basestring) and levs == 'all':
-            levs = numpy.array(range(maxLev+1))
-        elif isinstance(levs, basestring):
-            print "Error: %s not valid for levs parameter." % (levs)
-            print "levs must be either a 1D numpy array or the string 'all'."
-            return
+        if isinstance(levs, str) and levs == 'all':
+            levs = numpy.array(list(range(maxLev+1)))
+        elif isinstance(levs, str):
+            raise ValueError("Error: %s not valid for levs parameter. "
+                             "levs must be either a 1D numpy array or the string 'all'." % levs)
         else:
             levs = numpy.array(levs)
 
-        if isinstance(bands, basestring) and bands == 'all':
-            bands = numpy.array(range(nbands))
-        elif isinstance(bands, basestring):
-            print "Error: %s not valid for bands parameter." % (bands)
-            print "bands must be either a 1D numpy array or the string 'all'."
-            return
+        if isinstance(bands, str) and bands == 'all':
+            bands = numpy.array(list(range(nbands)))
+        elif isinstance(bands, str):
+            raise ValueError("Error: %s not valid for bands parameter. " \
+                  "bands must be either a 1D numpy array or the string 'all'." % bands)
         else:
             bands = numpy.array(bands)
 
@@ -224,30 +203,30 @@ def reconSFpyr(self, *args):
             dims = numpy.array(self.pyrSize[dimIdx])
             if (dims[0], dims[1]) not in dimList:
                 dimList.append( (dims[0], dims[1]) )
-            ctr = numpy.ceil((dims+0.5)/2)
-            lodims = numpy.ceil((dims-0.5)/2)
-            loctr = numpy.ceil((lodims+0.5)/2)
+            ctr = numpy.ceil((dims+0.5)/2).astype('int')
+            lodims = numpy.ceil((dims-0.5)/2).astype('int')
+            loctr = numpy.ceil((lodims+0.5)/2).astype('int')
             lostart = ctr - loctr
             loend = lostart + lodims
             bounds = (lostart[0], lostart[1], loend[0], loend[1])
             if bounds not in boundList:
                 boundList.append( bounds )
-        boundList.append((0.0, 0.0, dimList[len(dimList)-1][0], 
+        boundList.append((0, 0, dimList[len(dimList)-1][0],
                           dimList[len(dimList)-1][1]))
         dimList.append((dimList[len(dimList)-1][0], dimList[len(dimList)-1][1]))
         
         # matlab code starts here
         dims = numpy.array(self.pyrSize[0])
-        ctr = numpy.ceil((dims+0.5)/2.0)
+        ctr = numpy.ceil((dims+0.5)/2.0).astype('int')
 
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
-                                     (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
-                                     (dims[0]/2))
+        (xramp, yramp) = numpy.meshgrid((numpy.array(list(range(1,dims[1]+1)))-ctr[1])/
+                                     (dims[1]/2.),
+                                     (numpy.array(list(range(1,dims[0]+1)))-ctr[0])/
+                                     (dims[0]/2.))
         angle = numpy.arctan2(yramp, xramp)
         log_rad = numpy.sqrt(xramp**2 + yramp**2)
         log_rad[ctr[0]-1, ctr[1]-1] = log_rad[ctr[0]-1, ctr[1]-2]
-        log_rad = numpy.log2(log_rad);
+        log_rad = numpy.log2(log_rad)
 
         ## Radial transition function (a raised cosine in log-frequency):
         (Xrcos, Yrcos) = rcosFn(twidth, (-twidth/2.0), numpy.array([0,1]))
@@ -256,7 +235,7 @@ def reconSFpyr(self, *args):
 
         # from reconSFpyrLevs
         lutsize = 1024
-        Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
+        Xcosn = numpy.pi * numpy.array(list(range(-(2*lutsize+1), (lutsize+2)))) / lutsize
         
         order = nbands -1
         const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
@@ -278,6 +257,7 @@ def reconSFpyr(self, *args):
                       bounds[0]+boundList[idx][0] + diff[0], 
                       bounds[1]+boundList[idx][1] + diff[1])
             Xrcos -= numpy.log2(2.0)
+        bounds = numpy.array(bounds, dtype='int')
         nlog_rad = log_rad[bounds[0]:bounds[2], bounds[1]:bounds[3]]
 
         nlog_rad_tmp = numpy.reshape(nlog_rad, 
diff --git a/Spyr.py b/Spyr.py
index 554e39a..b3beee3 100644
--- a/Spyr.py
+++ b/Spyr.py
@@ -1,18 +1,18 @@
-from pyramid import pyramid
+from .pyramid import pyramid
 import numpy
-from sp0Filters import sp0Filters
-from sp1Filters import sp1Filters
-from sp3Filters import sp3Filters
-from sp5Filters import sp5Filters
+from .sp0Filters import sp0Filters
+from .sp1Filters import sp1Filters
+from .sp3Filters import sp3Filters
+from .sp5Filters import sp5Filters
 import os
-from maxPyrHt import maxPyrHt
-from corrDn import corrDn
+from .maxPyrHt import maxPyrHt
+from .corrDn import corrDn
 import math
-from LB2idx import LB2idx
+from .LB2idx import LB2idx
 import matplotlib
-from showIm import showIm
-import JBhelpers
-from upConv import upConv
+from .showIm import showIm
+from . import JBhelpers
+from .upConv import upConv
 
 class Spyr(pyramid):
     filt = ''
@@ -24,7 +24,7 @@ def __init__(self, *args):    # (image height, filter file, edges)
         if len(args) > 0:
             self.image = numpy.array(args[0])
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -40,10 +40,10 @@ def __init__(self, *args):    # (image height, filter file, edges)
             elif args[2] == 'sp5Filters':
                 filters = sp5Filters()
             elif os.path.isfile(args[2]):
-                print "Filter files not supported yet"
+                print("Filter files not supported yet")
                 return
             else:
-                print "filter parameters value %s not supported" % (args[2])
+                print("filter parameters value %s not supported" % (args[2]))
                 return
         else:
             filters = sp1Filters()
@@ -60,8 +60,8 @@ def __init__(self, *args):    # (image height, filter file, edges)
             if args[1] == 'auto':
                 ht = max_ht
             elif args[1] > max_ht:
-                print "Error: cannot build pyramid higher than %d levels." % (
-                    max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (
+                    max_ht))
                 return
             else:
                 ht = args[1]
@@ -115,20 +115,20 @@ def __init__(self, *args):    # (image height, filter file, edges)
     # methods
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        if isinstance(args[1], (int, long)):
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        if isinstance(args[1], int):
             self.pyr[args[0]][0][args[1]] = args[2]
         elif isinstance(args[1], tuple):
             self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
         else:
-            print 'Error: location parameter must be int or tuple!'
+            print('Error: location parameter must be int or tuple!')
             return
 
     def spyrLev(self, lev):
         if lev < 0 or lev > self.spyrHt()-1:
-            print 'Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1)
+            print('Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1))
             return
         
         levArray = []
@@ -140,13 +140,21 @@ def spyrLev(self, lev):
 
     def spyrBand(self, lev, band):
         if lev < 0 or lev > self.spyrHt()-1:
-            print 'Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1)
+            print('Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1))
             return
         if band < 0 or band > self.numBands()-1:
-            print 'Error: band parameter must be between 0 and %d!' % (self.numBands()-1)
+            print('Error: band parameter must be between 0 and %d!' % (self.numBands()-1))
 
         return self.band( ((lev*self.numBands())+band)+1 )
 
+    def bandIndex(self, lev, band):
+        if lev < 0 or lev > self.spyrHt()-1:
+            raise ValueError('Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1))
+        if band < 0 or band > self.numBands()-1:
+            raise ValueError('Error: band parameter must be between 0 and %d!' % (self.numBands()-1))
+
+        return ((lev*self.numBands())+band)+1
+
     def spyrHt(self):
         if len(self.pyrSize) > 2:
             spHt = (len(self.pyrSize)-2)/self.numBands()
@@ -183,10 +191,10 @@ def reconPyr(self, *args):
             elif args[0] == 'sp5Filters':
                 filters = sp5Filters()
             elif os.path.isfile(args[0]):
-                print "Filter files not supported yet"
+                print("Filter files not supported yet")
                 return
             else:
-                print "filter %s not supported" % (args[0])
+                print("filter %s not supported" % (args[0]))
                 return
         else:
             filters = sp1Filters()
@@ -218,22 +226,22 @@ def reconPyr(self, *args):
         
         maxLev = 2 + self.spyrHt()
         if levs == 'all':
-            levs = numpy.array(range(maxLev))
+            levs = numpy.array(list(range(maxLev)))
         else:
             levs = numpy.array(levs)
             if (levs < 0).any() or (levs >= maxLev).any():
-                print "Error: level numbers must be in the range [0, %d]." % (maxLev-1)
+                print("Error: level numbers must be in the range [0, %d]." % (maxLev-1))
                 return
             else:
                 levs = numpy.array(levs)
                 if len(levs) > 1 and levs[0] < levs[1]:
                     levs = levs[::-1]  # we want smallest first
         if bands == 'all':
-            bands = numpy.array(range(self.numBands()))
+            bands = numpy.array(list(range(self.numBands())))
         else:
             bands = numpy.array(bands)
             if (bands < 0).any() or (bands > bfilts.shape[1]).any():
-                print "Error: band numbers must be in the range [0, %d]." % (self.numBands()-1)
+                print("Error: band numbers must be in the range [0, %d]." % (self.numBands()-1))
                 return
             else:
                 bands = numpy.array(bands)
@@ -359,10 +367,10 @@ def showPyr(self, prange = 'auto2', gap = 1, scale = 2, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.sqrt( numpy.var(band) )
             prange[nind-1,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(prange, basestring):
-            print "Error:Bad RANGE argument: %s'" % (prange)
+        elif isinstance(prange, str):
+            print("Error:Bad RANGE argument: %s'" % (prange))
         elif prange.shape[0] == 1 and prange.shape[1] == 2:
-            scales = numpy.power(scale, range(ht))
+            scales = numpy.power(scale, list(range(ht)))
             scales = numpy.outer( numpy.ones((nbands,1)), scales )
             scales = numpy.array([1, scales, numpy.power(scale, ht)])
             prange = numpy.outer(scales, prange)
@@ -381,11 +389,11 @@ def showPyr(self, prange = 'auto2', gap = 1, scale = 2, disp = 'qt'):
             ncols = int(numpy.ceil((nbands+1)/2))
             nrows = int(numpy.ceil(nbands/2))
 
-        a = numpy.array(range(1-nrows, 1))
+        a = numpy.array(list(range(1-nrows, 1)))
         b = numpy.zeros((1,ncols))[0]
         ab = numpy.concatenate((a,b))
         c = numpy.zeros((1,nrows))[0]
-        d = range(-1, -ncols-1, -1)
+        d = list(range(-1, -ncols-1, -1))
         cd = numpy.concatenate((c,d))
         relpos = numpy.vstack((ab,cd)).T
         
diff --git a/TESTS/unitTests.py b/TESTS/unitTests.py
index aaea27b..f297a0d 100755
--- a/TESTS/unitTests.py
+++ b/TESTS/unitTests.py
@@ -11,41 +11,41 @@
 
 class maxPyrHtTests(unittest.TestCase):
     def test1(self):
-        self.failUnless(ppt.maxPyrHt((1,10),(3,4)) == 0)
+        self.assertTrue(ppt.maxPyrHt((1,10),(3,4)) == 0)
     def test2(self):
-        self.failUnless(ppt.maxPyrHt((10,1),(3,4)) == 0)
+        self.assertTrue(ppt.maxPyrHt((10,1),(3,4)) == 0)
     def test3(self):
-        self.failUnless(ppt.maxPyrHt((10,10),(1,4)) == 4)
+        self.assertTrue(ppt.maxPyrHt((10,10),(1,4)) == 4)
     def test4(self):
-        self.failUnless(ppt.maxPyrHt((10,10),(3,1)) == 2)
+        self.assertTrue(ppt.maxPyrHt((10,10),(3,1)) == 2)
     def test5(self):
-        self.failUnless(ppt.maxPyrHt((10,10),(3,4)) == 2)
+        self.assertTrue(ppt.maxPyrHt((10,10),(3,4)) == 2)
     def test6(self):
-        self.failUnless(ppt.maxPyrHt((20,10),(5,1)) == 2)
+        self.assertTrue(ppt.maxPyrHt((20,10),(5,1)) == 2)
     def test7(self):
-        self.failUnless(ppt.maxPyrHt((10,20),(5,1)) == 2)
+        self.assertTrue(ppt.maxPyrHt((10,20),(5,1)) == 2)
     def test8(self):
-        self.failUnless(ppt.maxPyrHt((20,10),(1,5)) == 5)
+        self.assertTrue(ppt.maxPyrHt((20,10),(1,5)) == 5)
     def test9(self):
-        self.failUnless(ppt.maxPyrHt((10,20),(1,5)) == 5)
+        self.assertTrue(ppt.maxPyrHt((10,20),(1,5)) == 5)
     def test10(self):
-        self.failUnless(ppt.maxPyrHt((256,1),(1,5)) == 6)
+        self.assertTrue(ppt.maxPyrHt((256,1),(1,5)) == 6)
     def test11(self):
-        self.failUnless(ppt.maxPyrHt((256,1),(5,1)) == 6)
+        self.assertTrue(ppt.maxPyrHt((256,1),(5,1)) == 6)
     def test12(self):
-        self.failUnless(ppt.maxPyrHt((1,256),(1,5)) == 6)
+        self.assertTrue(ppt.maxPyrHt((1,256),(1,5)) == 6)
     def test13(self):
-        self.failUnless(ppt.maxPyrHt((1,256),(5,1)) == 6)
+        self.assertTrue(ppt.maxPyrHt((1,256),(5,1)) == 6)
 
 class binomialFilterTests(unittest.TestCase):
     def test1(self):
-        self.failUnless((ppt.binomialFilter(2) == np.array([[0.5],
+        self.assertTrue((ppt.binomialFilter(2) == np.array([[0.5],
                                                             [0.5]])).all() )
     def test2(self):
-        self.failUnless((ppt.binomialFilter(3) == np.array([[0.25], [0.5],
+        self.assertTrue((ppt.binomialFilter(3) == np.array([[0.25], [0.5],
                                                             [0.25]])).all())
     def test3(self):
-        self.failUnless((ppt.binomialFilter(5) == np.array([[0.0625], [0.25],
+        self.assertTrue((ppt.binomialFilter(5) == np.array([[0.0625], [0.25],
                                                             [0.3750], [0.25],
                                                             [0.0625]])).all())
 
@@ -55,34 +55,34 @@ def test1(self):
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/buildGpyr2row.mat')
-        img = np.array(range(256)).astype(float)
+        img = np.array(list(range(256))).astype(float)
         img = img.reshape(1, 256)
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/buildGpyr2col.mat')
-        img = np.array(range(256)).astype(float)
+        img = np.array(list(range(256))).astype(float)
         img = img.reshape(256, 1)
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildGpyr3.mat')
         img = ppt.mkRamp(10)
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/buildGpyr4.mat')
         img = ppt.mkRamp((10,20))
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildGpyr5.mat')
         img = ppt.mkRamp((20, 10))
         pyPyr = ppt.Gpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
 
 class LpyrTests(unittest.TestCase):
     def test1(self):
@@ -90,289 +90,289 @@ def test1(self):
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.Lpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self): 
         matPyr2 = scipy.io.loadmat('../matFiles/buildLpyr2.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Lpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr2['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr2['pyr'], pyPyr))
     def test3(self):
         matPyr2 = scipy.io.loadmat('../matFiles/buildLpyr3.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Lpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr2['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr2['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr4.mat')
         pyRamp = ppt.mkRamp(200,100)
         pyPyr = ppt.Lpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test5(self):   
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr5.mat')
-        pyRamp = np.array(range(200)).reshape(1, 200)
+        pyRamp = np.array(list(range(200))).reshape(1, 200)
         pyPyr = ppt.Lpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr6.mat')
-        pyRamp = np.array(range(200))
+        pyRamp = np.array(list(range(200)))
         pyPyr = ppt.Lpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr7.mat')
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.Lpyr(img)
         recon = pyPyr.reconPyr()
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
     def test8(self): 
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr8.mat')
         pyRamp = ppt.mkRamp(200)
         pyPyr = ppt.Lpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
     def test9(self): 
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr9.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.Lpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
     def test10(self): 
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr10.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Lpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
     def test11(self): 
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr11.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Lpyr(pyRamp)
         recon = pyPyr.reconPyr([1])
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
     def test12(self): 
         matPyr = scipy.io.loadmat('../matFiles/buildLpyr12.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Lpyr(pyRamp)
         recon = pyPyr.reconPyr([0,2,4])
-        self.failUnless((matPyr['recon'] == recon).all())
+        self.assertTrue((matPyr['recon'] == recon).all())
 
 class spFilterTests(unittest.TestCase):
     def test1(self):
         matFilt0 = scipy.io.loadmat('../matFiles/sp0Filters.mat')
         pySP0filt = ppt.sp0Filters()
         tmpKeys = []
-        for key in matFilt0.keys():
+        for key in list(matFilt0.keys()):
             if "_" not in key:
                 tmpKeys.append(key)
-        self.failUnless(tmpKeys == pySP0filt.keys())
+        self.assertTrue(tmpKeys == list(pySP0filt.keys()))
         for key in tmpKeys:
-            self.failUnless((matFilt0[key] == pySP0filt[key]).all())
+            self.assertTrue((matFilt0[key] == pySP0filt[key]).all())
 
     def test2(self):
         matFilt1 = scipy.io.loadmat('../matFiles/sp1Filters.mat')
         pySP1filt = ppt.sp1Filters()
         tmpKeys = []
-        for key in matFilt1.keys():
+        for key in list(matFilt1.keys()):
             if "_" not in key:
                 tmpKeys.append(key)
-        self.failUnless(tmpKeys == pySP1filt.keys())
+        self.assertTrue(tmpKeys == list(pySP1filt.keys()))
         for key in tmpKeys:
-            self.failUnless((matFilt1[key] == pySP1filt[key]).all())
+            self.assertTrue((matFilt1[key] == pySP1filt[key]).all())
 
     def test3(self):
         matFilt3 = scipy.io.loadmat('../matFiles/sp3Filters.mat')
         pySP3filt = ppt.sp3Filters()
         tmpKeys = []
-        for key in matFilt3.keys():
+        for key in list(matFilt3.keys()):
             if "_" not in key:
                 tmpKeys.append(key)
-        self.failUnless(tmpKeys == pySP3filt.keys())
+        self.assertTrue(tmpKeys == list(pySP3filt.keys()))
         for key in tmpKeys:
-            self.failUnless((matFilt3[key] == pySP3filt[key]).all())
+            self.assertTrue((matFilt3[key] == pySP3filt[key]).all())
 
     def test4(self):
         matFilt5 = scipy.io.loadmat('../matFiles/sp5Filters.mat')
         pySP5filt = ppt.sp5Filters()
         tmpKeys = []
-        for key in matFilt5.keys():
+        for key in list(matFilt5.keys()):
             if "_" not in key:
                 tmpKeys.append(key)
-        self.failUnless(tmpKeys == pySP5filt.keys())
+        self.assertTrue(tmpKeys == list(pySP5filt.keys()))
         for key in tmpKeys:
-            self.failUnless((matFilt5[key] == pySP5filt[key]).all())
+            self.assertTrue((matFilt5[key] == pySP5filt[key]).all())
 
 class SpyrTests(unittest.TestCase):
     def test00(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr00.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.Spyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr1.mat')
         img = np.array(Image.open('../lenna-256x256.tif')).astype(float)
         pyPyr = ppt.Spyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr2.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr3.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Spyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr4.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.Spyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr0.mat')
         pyRamp = ppt.mkRamp(20)
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr5.mat')
         img = np.array(Image.open('../lenna-256x256.tif')).astype(float)
         pyPyr = ppt.Spyr(img)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr6.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr7.mat')
         pyRamp = ppt.mkRamp((256,128))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test8(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr8.mat')
         pyRamp = ppt.mkRamp((128,256))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test9(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr9.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test10(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr10.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test11(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr11.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr('sp1Filters', 'reflect1', [1])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test12(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr12.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp)
         recon = pyPyr.reconPyr('sp1Filters', 'reflect1', [0,2,4])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test13(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr13.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.Spyr(pyRamp, 1, 'sp0Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test14(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr14.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 3, 'sp0Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test15(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr15.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 1, 'sp1Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test16(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr16.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 3, 'sp1Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test17(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr17.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 1, 'sp3Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test18(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr18.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 3, 'sp3Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test19(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr19.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 1, 'sp5Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test20(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr20.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Spyr(pyRamp, 3, 'sp5Filters')
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test21(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr21.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:256,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp0Filters')
         recon = pyPyr.reconPyr('sp0Filters');
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test22(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr22.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:256,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp1Filters')
         recon = pyPyr.reconPyr('sp1Filters');
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test23(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr23.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:256,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp3Filters')
         recon = pyPyr.reconPyr('sp3Filters');
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test24(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr24.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:256,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp5Filters')
         recon = pyPyr.reconPyr('sp5Filters');
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test25(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr25.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:256,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp5Filters')
         recon = pyPyr.reconPyr('sp5Filters','reflect1',[0,1,2], [0]);
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test26(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr26.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:128,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp0Filters')
         recon = pyPyr.reconPyr('sp0Filters','reflect1',[0,1,2], [0]);
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test27(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr27.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:128,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp1Filters')
         recon = pyPyr.reconPyr('sp1Filters','reflect1',[0,1,2], [0]);
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test28(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr28.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:128,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp3Filters')
         recon = pyPyr.reconPyr('sp3Filters','reflect1',[0,1,2], [0]);
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test29(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSpyr29.mat')
         texture = scipy.io.loadmat('../matFiles/im04-1.mat')['res'][0:128,0:256]
         pyPyr = ppt.Spyr(texture, 3, 'sp5Filters')
         recon = pyPyr.reconPyr('sp5Filters','reflect1',[0,1,2], [0]);
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
         
         
         
@@ -385,137 +385,137 @@ def test1(self):
         #foo = pointOp(200, 200, img, 5, filt, 0, 1, 0);
         foo = ppt.pointOp(img, filt, 0, 1, 0);
         foo = np.reshape(foo,(200,200))
-        self.failUnless((matImg['foo'] == foo).all())
+        self.assertTrue((matImg['foo'] == foo).all())
 
 class SFpyrTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr0.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.SFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr1.mat')
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.SFpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr2.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.SFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr3.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.SFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr4.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.SFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr5.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr6.mat')
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.SFpyr(img)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr7.mat')
         pyRamp = ppt.mkRamp((256,128))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test8(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr8.mat')
         pyRamp = ppt.mkRamp((128,256))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test9(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr9.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test10(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr10.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test11(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr11.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr([0])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test12(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr12.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr([0,2,4])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test13(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSFpyr13.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.SFpyr(pyRamp)
         recon = pyPyr.reconPyr([0,2,4], [1])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
 
 class SCFpyrTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr0.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.SCFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr1.mat')
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.SCFpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr2.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.SCFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr3.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.SCFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr4.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.SCFpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr5.mat')
         img = Image.open('../lenna-256x256.tif')
         img = np.array(img.getdata()).reshape(256,256)
         pyPyr = ppt.SCFpyr(img)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr6.mat')
         pyRamp = ppt.mkRamp((256,128))
         pyPyr = ppt.SCFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr7.mat')
         pyRamp = ppt.mkRamp((128,256))
         pyPyr = ppt.SCFpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     #def test8(self):  # fails in matlab version
     #    matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr8.mat')
     #    pyRamp = ppt.mkRamp((200,100))
@@ -533,170 +533,170 @@ def test10(self):
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.SCFpyr(pyRamp)
         recon = pyPyr.reconPyr([0])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test11(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr11.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.SCFpyr(pyRamp)
         recon = pyPyr.reconPyr([0,2,4])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test12(self):
         matPyr = scipy.io.loadmat('../matFiles/buildSCFpyr12.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.SCFpyr(pyRamp)
         recon = pyPyr.reconPyr([0,2,4], [1])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
 
 class WpyrTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr0.mat')
         pyRamp = ppt.mkRamp((20,20))
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr1.mat')
         img = np.array(Image.open('../lenna-256x256.tif')).astype(float)
         pyPyr = ppt.Wpyr(img)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr2.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr3.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr4.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr5.mat')
         img = np.array(Image.open('../lenna-256x256.tif')).astype(float)
         pyPyr = ppt.Wpyr(img)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr6.mat')
         pyRamp = ppt.mkRamp((256,128))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr7.mat')
         pyRamp = ppt.mkRamp((128,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test8(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr8.mat')
         pyRamp = ppt.mkRamp((200,200))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test9(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr9.mat')
         pyRamp = ppt.mkRamp((200,100))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test10(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr10.mat')
         pyRamp = ppt.mkRamp((100,200))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr()
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test11(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr11.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf9', 'reflect1', [0])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test12(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr12.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf9', 'reflect1', [0,2,4])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test13(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr13.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf9', 'reflect1', [0,2,4], [1])
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test14(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr14.mat')
         pyRamp = ppt.mkRamp((256,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf8')
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test15(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr15.mat')
         pyRamp = ppt.mkRamp((256,128))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf8')
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test16(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr16.mat')
         pyRamp = ppt.mkRamp((128,256))
         pyPyr = ppt.Wpyr(pyRamp)
         recon = pyPyr.reconPyr('qmf8')
-        self.failUnless(ppt.compareRecon(matPyr['recon'], recon))
+        self.assertTrue(ppt.compareRecon(matPyr['recon'], recon))
     def test17(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr17.mat')
         pyRamp = ppt.mkRamp((1,200))
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
     def test18(self):
         matPyr = scipy.io.loadmat('../matFiles/buildWpyr18.mat')
         pyRamp = ppt.mkRamp((1,200)).T
         pyPyr = ppt.Wpyr(pyRamp)
-        self.failUnless(ppt.comparePyr(matPyr['pyr'], pyPyr))
+        self.assertTrue(ppt.comparePyr(matPyr['pyr'], pyPyr))
 
 class blurDnTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn0.mat')
         pyRamp = ppt.mkRamp((20,20))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn1.mat')
         pyRamp = ppt.mkRamp((256,256))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn2.mat')
         pyRamp = ppt.mkRamp((256,128))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn3.mat')
         pyRamp = ppt.mkRamp((128,256))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn4.mat')
         pyRamp = ppt.mkRamp((200, 100))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn5.mat')
         pyRamp = ppt.mkRamp((100, 200))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn6.mat')
         pyRamp = ppt.mkRamp((1, 256)).T
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/blurDn7.mat')
         pyRamp = ppt.mkRamp((1, 256))
         res = ppt.blurDn(pyRamp)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     #def test8(self):  need a 2D filter
     #    matPyr = scipy.io.loadmat('../matFiles/blurDn8.mat')
     #    pyRamp = ppt.mkRamp((256, 256))
@@ -707,243 +707,243 @@ class mkAngularSineTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkAngularSine0.mat')
         res = ppt.mkAngularSine(20)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkAngularSine1.mat')
         res = ppt.mkAngularSine(20, 5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkAngularSine2.mat')
         res = ppt.mkAngularSine(20, 5, 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/mkAngularSine3.mat')
         res = ppt.mkAngularSine(20, 5, 3, 2)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/mkAngularSine4.mat')
         res = ppt.mkAngularSine(20, 5, 3, 2, (2,2))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class mkGaussianTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkGaussian0.mat')
         res = ppt.mkGaussian(20)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkGaussian1.mat')
         res = ppt.mkGaussian(20, (2,3))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkGaussian2.mat')
         res = ppt.mkGaussian(20, [[-1, 0], [0, 1]])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/mkGaussian3.mat')
         res = ppt.mkGaussian(10, [[-1, 0], [0, 1]])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/mkGaussian4.mat')
         res = ppt.mkGaussian(20, [[2, 0], [0, 1]])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class mkDiscTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkDisc0.mat')
         res = ppt.mkDisc(20)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkDisc1.mat')
         res = ppt.mkDisc(20, 8)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkDisc2.mat')
         res = ppt.mkDisc(20, 8, (0,0))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/mkDisc3.mat')
         res = ppt.mkDisc(20, 8, (0,0), 5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/mkDisc4.mat')
         res = ppt.mkDisc(20, 8, (0,0), 5, (0.75, 0.25))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class mkSineTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine0.mat')
         res = ppt.mkSine(20, 5.5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine1.mat')
         res = ppt.mkSine(20, 5.5, 2)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine2.mat')
         res = ppt.mkSine(20, 5.5, 2, 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine3.mat')
         res = ppt.mkSine(20, 5.5, 2, 3, 5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine4.mat')
         res = ppt.mkSine(20, 5.5, 2, 3, 5, [4,5])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine5.mat')
         res = ppt.mkSine(20, [1,2])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine6.mat')
         res = ppt.mkSine(20, [1,2], 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine7.mat')
         res = ppt.mkSine(20, [1,2], 3, 2)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test8(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSine8.mat')
         res = ppt.mkSine(20, [1,2], 3, 2, [5,4])
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class mkZonePlateTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkZonePlate0.mat')
         res = ppt.mkZonePlate(20)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkZonePlate1.mat')
         res = ppt.mkZonePlate(20, 4)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkZonePlate2.mat')
         res = ppt.mkZonePlate(20, 4, 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class mkSquareTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare0.mat')
         res = ppt.mkSquare(20, 5.5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare1.mat')
         res = ppt.mkSquare(20, 5.5, 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare2.mat')
         res = ppt.mkSquare(20, 5.5, 3, 5.1)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare3.mat')
         res = ppt.mkSquare(20, 5.5, 3, 5.1, -1)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare4.mat')
         res = ppt.mkSquare(20, 5.5, 3, 5.1, -1, (2,3))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare5.mat')
         res = ppt.mkSquare(20, 5.5, 3, 5.1, -1, (2,3), 0.25)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare6.mat')
         res = ppt.mkSquare(20, (1,2))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare7.mat')
         res = ppt.mkSquare(20, (1,2), 3.2)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test8(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare8.mat')
         res = ppt.mkSquare(20, (1,2), 3.2, -2)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test9(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare9.mat')
         res = ppt.mkSquare(20, (1,2), 3.2, -2, (2,3))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test10(self):
         matPyr = scipy.io.loadmat('../matFiles/mkSquare10.mat')
         res = ppt.mkSquare(20, (1,2), 3.2, -2, (2,3), 0.55)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class blurTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/blur0.mat')
         res = ppt.blur(ppt.mkRamp(20))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/blur1.mat')
         res = ppt.blur(ppt.mkRamp(20), 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/blur2.mat')
         res = ppt.blur(ppt.mkRamp(20), 3, ppt.namedFilter('qmf5'))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/blur3.mat')
         res = ppt.blur(ppt.mkRamp((20,30)))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/blur4.mat')
         res = ppt.blur(ppt.mkRamp((20,30)), 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/blur5.mat')
         res = ppt.blur(ppt.mkRamp((20,30)), 3, ppt.namedFilter('qmf5'))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class cconv2Tests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_0.mat')
         res = ppt.cconv2(ppt.mkRamp(20), ppt.mkRamp(10))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_1.mat')
         res = ppt.cconv2(ppt.mkRamp(10), ppt.mkRamp(20))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_2.mat')
         res = ppt.cconv2(ppt.mkRamp(20), ppt.mkRamp(10), 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_3.mat')
         res = ppt.cconv2(ppt.mkRamp(10), ppt.mkRamp(20), 3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_4.mat')
         res = ppt.cconv2(ppt.mkRamp((20,30)), ppt.mkRamp((10,20)))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_5.mat')
         res = ppt.cconv2(ppt.mkRamp((10,20)), ppt.mkRamp((20,30)))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_6.mat')
         res = ppt.cconv2(ppt.mkRamp((20,30)), ppt.mkRamp((10,20)), 5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test7(self):
         matPyr = scipy.io.loadmat('../matFiles/cconv2_7.mat')
         res = ppt.cconv2(ppt.mkRamp((10,20)), ppt.mkRamp((20,30)), 5)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 class clipTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/clip0.mat')
         res = ppt.clip(ppt.mkRamp(20) / ppt.mkRamp(20).max())
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/clip1.mat')
         res = ppt.clip(ppt.mkRamp(20) / ppt.mkRamp(20).max(), 0.3)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/clip2.mat')
         res = ppt.clip(ppt.mkRamp(20) / ppt.mkRamp(20).max(), (0.3,0.7))
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/clip3.mat')
         res = ppt.clip(ppt.mkRamp(20) / ppt.mkRamp(20).max(), 0.3, 0.7)
-        self.failUnless(ppt.compareRecon(matPyr['res'], res))
+        self.assertTrue(ppt.compareRecon(matPyr['res'], res))
 
 # python version of histo
 # adding 0.7 to ramp to nullify rounding differences between Python and Matlab
@@ -955,7 +955,7 @@ def test0(self):
         #   python return edges
         #self.failUnless(ppt.compareRecon(matPyr['X'], X))
         #self.failUnless(ppt.compareRecon(matPyr['N'], N))
-        self.failUnless((matPyr['N'] == N).all())
+        self.assertTrue((matPyr['N'] == N).all())
 # FIX: why does matlab version return N+1 bins??
 #    def test1(self):
 #        matPyr = scipy.io.loadmat('../matFiles/histo1.mat')
@@ -994,11 +994,11 @@ class factorialTests(unittest.TestCase):
     def test0(self):
         matPyr = scipy.io.loadmat('../matFiles/factorial0.mat')
         res = ppt.factorial([[1,2],[3,4]])
-        self.failUnless((matPyr['res'] == res).all())
+        self.assertTrue((matPyr['res'] == res).all())
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/factorial1.mat')
         res = ppt.factorial(4)
-        self.failUnless(matPyr['res'] == res)
+        self.assertTrue(matPyr['res'] == res)
 
 #class histoMatchTests(unittest.TestCase):
 #    def test0(self):
@@ -1020,16 +1020,16 @@ def test0(self):
         [dx,dy] = ppt.imGradient(ramp)
         dx = np.array(dx)
         dy = np.array(dy)
-        self.failUnless(ppt.compareRecon(matPyr['res'][:,:,0], dx))
-        self.failUnless(ppt.compareRecon(matPyr['res'][:,:,1], dy))
+        self.assertTrue(ppt.compareRecon(matPyr['res'][:,:,0], dx))
+        self.assertTrue(ppt.compareRecon(matPyr['res'][:,:,1], dy))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/imGradient1.mat')
         ramp = ppt.mkRamp(10)
         [dx,dy] = ppt.imGradient(ramp, 'reflect1')
         dx = np.array(dx)
         dy = np.array(dy)
-        self.failUnless(ppt.compareRecon(matPyr['res'][:,:,0], dx))
-        self.failUnless(ppt.compareRecon(matPyr['res'][:,:,1], dy))
+        self.assertTrue(ppt.compareRecon(matPyr['res'][:,:,0], dx))
+        self.assertTrue(ppt.compareRecon(matPyr['res'][:,:,1], dy))
  
 class skew2Tests(unittest.TestCase):
     def test0(self):
@@ -1038,7 +1038,7 @@ def test0(self):
         mres = matPyr['res'][0][0]
         disc = ppt.mkDisc(10)
         res = ppt.skew2(disc)
-        self.failUnless(np.absolute(res - mres) <= math.pow(10,-11))
+        self.assertTrue(np.absolute(res - mres) <= math.pow(10,-11))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/skew2_1.mat')
         # not sure why matPyr is [[ans]]???
@@ -1047,7 +1047,7 @@ def test1(self):
         # using incorrect mean for better test
         mn = disc.mean() + 0.1
         res = ppt.skew2(disc, mn)
-        self.failUnless(np.absolute(res - mres) <= math.pow(10,-11))
+        self.assertTrue(np.absolute(res - mres) <= math.pow(10,-11))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/skew2_2.mat')
         # not sure why matPyr is [[ans]]???
@@ -1057,7 +1057,7 @@ def test2(self):
         mn = disc.mean() + 0.1
         v = ppt.var2(disc) + 0.1
         res = ppt.skew2(disc, mn, v)
-        self.failUnless(np.absolute(res - mres) <= math.pow(10,-11))
+        self.assertTrue(np.absolute(res - mres) <= math.pow(10,-11))
        
 class upBlurTests(unittest.TestCase):
     def test0(self):
@@ -1065,44 +1065,44 @@ def test0(self):
         mres = matPyr['res']
         im = ppt.mkRamp((1,20))
         res = ppt.upBlur(im)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur1.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((1,20))
         res = ppt.upBlur(im.T)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur2.mat')
         mres = matPyr['res']
         im = ppt.mkRamp(20)
         res = ppt.upBlur(im)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur3.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((1,20))
         res = ppt.upBlur(im, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur4.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((1,20))
         res = ppt.upBlur(im.T, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur5.mat')
         mres = matPyr['res']
         im = ppt.mkRamp(20)
         res = ppt.upBlur(im, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test6(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur6.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((1,20))
         filt = ppt.namedFilter('qmf9')
         res = ppt.upBlur(im, 3, filt)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     #def test7(self):   # fails in matlab and python because of dim mismatch
     #    matPyr = scipy.io.loadmat('../matFiles/upBlur7.mat')
     #    mres = matPyr['res']
@@ -1116,7 +1116,7 @@ def test8(self):
         im = ppt.mkRamp((1,20))
         filt = ppt.namedFilter('qmf9')
         res = ppt.upBlur(im.T, 3, filt)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     #def test9(self):  # fails in matlab and python because of dim mismatch
     #    matPyr = scipy.io.loadmat('../matFiles/upBlur6.mat')
     #    mres = matPyr['res']
@@ -1130,21 +1130,21 @@ def test10(self):
         im = ppt.mkRamp(20)
         filt = ppt.mkDisc(3)
         res = ppt.upBlur(im, 3, filt)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test11(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur11.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((20,10))
         filt = ppt.mkDisc((5,3))
         res = ppt.upBlur(im, 3, filt)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test12(self):
         matPyr = scipy.io.loadmat('../matFiles/upBlur12.mat')
         mres = matPyr['res']
         im = ppt.mkRamp((10,20))
         filt = ppt.mkDisc((3,5))
         res = ppt.upBlur(im, 3, filt)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     
 class zconv2Tests(unittest.TestCase):
     def test0(self):
@@ -1153,42 +1153,42 @@ def test0(self):
         ramp = ppt.mkRamp(10)
         disc = ppt.mkDisc(5)
         res = ppt.zconv2(ramp, disc)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test1(self):
         matPyr = scipy.io.loadmat('../matFiles/zconv2_1.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp((10,20))
         disc = ppt.mkDisc((5,10))
         res = ppt.zconv2(ramp, disc)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/zconv2_2.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp((20,10))
         disc = ppt.mkDisc((10,5))
         res = ppt.zconv2(ramp, disc)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/zconv2_3.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp(10)
         disc = ppt.mkDisc(5)
         res = ppt.zconv2(ramp, disc, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test4(self):
         matPyr = scipy.io.loadmat('../matFiles/zconv2_4.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp((10,20))
         disc = ppt.mkDisc((5,10))
         res = ppt.zconv2(ramp, disc, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test5(self):
         matPyr = scipy.io.loadmat('../matFiles/zconv2_5.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp((20,10))
         disc = ppt.mkDisc((10,5))
         res = ppt.zconv2(ramp, disc, 3)
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
 
 class corrDnTests(unittest.TestCase):
     def test1(self):
@@ -1196,19 +1196,19 @@ def test1(self):
         mres = matPyr['res']
         ramp = ppt.mkRamp(20)
         res = ppt.corrDn(ramp, ppt.namedFilter('binom5'))
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test2(self):
         matPyr = scipy.io.loadmat('../matFiles/corrDn2.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp(20)
         res = ppt.corrDn(ramp, ppt.namedFilter('qmf13'))
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
     def test3(self):
         matPyr = scipy.io.loadmat('../matFiles/corrDn3.mat')
         mres = matPyr['res']
         ramp = ppt.mkRamp(20)
         res = ppt.corrDn(ramp, ppt.namedFilter('qmf16'))
-        self.failUnless(ppt.compareRecon(mres, res))
+        self.assertTrue(ppt.compareRecon(mres, res))
 
 def main():
     unittest.main()
diff --git a/Wpyr.py b/Wpyr.py
index 1b1f6ac..63acd1a 100644
--- a/Wpyr.py
+++ b/Wpyr.py
@@ -1,11 +1,11 @@
-from Lpyr import Lpyr
-from LB2idx import LB2idx
-from namedFilter import namedFilter
-from modulateFlip import modulateFlip
-from maxPyrHt import maxPyrHt
-from corrDn import corrDn
-from upConv import upConv
-import JBhelpers
+from .Lpyr import Lpyr
+from .LB2idx import LB2idx
+from .namedFilter import namedFilter
+from .modulateFlip import modulateFlip
+from .maxPyrHt import maxPyrHt
+from .corrDn import corrDn
+from .upConv import upConv
+from . import JBhelpers
 import numpy
 import matplotlib
 import pylab
@@ -24,7 +24,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         if len(args) > 0:
             im = args[0]
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -34,11 +34,11 @@ def __init__(self, *args):    # (image, height, order, twidth)
             filt = args[2]
         else:
             filt = "qmf9"
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = namedFilter(filt)
 
         if len(filt.shape) != 1 and filt.shape[0] != 1 and filt.shape[1] != 1:
-            print "Error: filter should be 1D (i.e., a vector)";
+            print("Error: filter should be 1D (i.e., a vector)");
             return
         hfilt = modulateFlip(filt)
 
@@ -67,7 +67,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             if ht == 'auto':
                 ht = max_ht
             elif(ht > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (max_ht))
         else:
             ht = max_ht
         ht = int(ht)
@@ -166,28 +166,28 @@ def reconPyr(self, *args):
         maxLev = int(self.wpyrHt() + 1)
 
         if isinstance(levs, str) and levs == 'all':
-            levs = numpy.array(range(maxLev))
+            levs = numpy.array(list(range(maxLev)))
         else:
             tmpLevs = []
             for l in levs:
                 tmpLevs.append((maxLev-1)-l)
             levs = numpy.array(tmpLevs)
             if (levs > maxLev).any():
-                print "Error: level numbers must be in the range [0, %d]" % (maxLev)
-        allLevs = numpy.array(range(maxLev))
+                print("Error: level numbers must be in the range [0, %d]" % (maxLev))
+        allLevs = numpy.array(list(range(maxLev)))
 
         if isinstance(bands, str) and bands == "all":
             if ( len(self.band(0)) == 1 or self.band(0).shape[0] == 1 or 
                  self.band(0).shape[1] == 1 ):
                 bands = numpy.array([0]);
             else:
-                bands = numpy.array(range(3))
+                bands = numpy.array(list(range(3)))
         else:
             bands = numpy.array(bands)
             if (bands < 0).any() or (bands > 2).any():
-                print "Error: band numbers must be in the range [0,2]."
+                print("Error: band numbers must be in the range [0,2].")
         
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = namedFilter(filt)
 
         hfilt = modulateFlip(filt).T
@@ -302,25 +302,25 @@ def reconPyr(self, *args):
 
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        if isinstance(args[1], (int, long)):
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        if isinstance(args[1], int):
             self.pyr[args[0]][0][args[1]] = args[2]
         elif isinstance(args[1], tuple):
             self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
         else:
-            print 'Error: location parameter must be int or tuple!'
+            print('Error: location parameter must be int or tuple!')
             return
             
 
     def set1D(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        print '%d %d %d' % (args[0], args[1], args[2])
-        print self.pyr[args[0]][0][1]
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        print('%d %d %d' % (args[0], args[1], args[2]))
+        print(self.pyr[args[0]][0][1])
 
     def pyrLow(self):
         return numpy.array(self.band(len(self.pyrSize)-1))
@@ -410,10 +410,10 @@ def showPyr(self, prange = None, gap = 1, scale = None, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.sqrt( numpy.var(band) )
             prange[nind-1,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(prange, basestring):
-            print "Error:Bad RANGE argument: %s'" % (prange)
+        elif isinstance(prange, str):
+            print("Error:Bad RANGE argument: %s'" % (prange))
         elif prange.shape[0] == 1 and prange.shape[1] == 2:
-            scales = numpy.power(scale, range(ht))
+            scales = numpy.power(scale, list(range(ht)))
             scales = numpy.outer( numpy.ones((nbands,1)), scales )
             scales = numpy.array([1, scales, numpy.power(scale, ht)])
             prange = numpy.outer(scales, prange)
diff --git a/__init__.py b/__init__.py
index 204161c..322cabf 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,66 +1,66 @@
-from binomialFilter import binomialFilter
-from blurDn import blurDn
-from blur import blur
-from cconv2 import cconv2
-from clip import clip
-from comparePyr import comparePyr
-from compareRecon import compareRecon
-from corrDn import corrDn
-from entropy2 import entropy2
-from factorial import factorial
-from Gpyr import Gpyr
-from histoMatch import histoMatch
-from histo import histo
-from idx2LB import idx2LB
-from imGradient import imGradient
-from imStats import imStats
-from kurt2 import kurt2
-from LB2idx import LB2idx
-from Lpyr import Lpyr
-from maxPyrHt import maxPyrHt
-from mkAngle import mkAngle
-from mkAngularSine import mkAngularSine
-from mkDisc import mkDisc
-from mkFract import mkFract
-from mkGaussian import mkGaussian
-from mkImpulse import mkImpulse
-from mkRamp import mkRamp
-from mkR import mkR
-from mkSine import mkSine
-from mkSquare import mkSquare
-from mkZonePlate import mkZonePlate
-from modulateFlip import modulateFlip
-from namedFilter import namedFilter
-from nextSz import nextSz
-from pointOp import pointOp
-from pyramid import pyramid
-from range2 import range2
-from rconv2 import rconv2
-from rcosFn import rcosFn
-from round import round
-from roundVal import roundVal
-from SCFpyr import SCFpyr
-from SFpyr import SFpyr
-from shift import shift
-from showIm import showIm
-from skew2 import skew2
-from sp0Filters import sp0Filters
-from sp1Filters import sp1Filters
-from sp3Filters import sp3Filters
-from sp5Filters import sp5Filters
-from Spyr import Spyr
-from SpyrMulti import SpyrMulti
-from steer2HarmMtx import steer2HarmMtx
-from steer import steer
-from strictly_decreasing import strictly_decreasing
-from upBlur import upBlur
-from upConv import upConv
-from var2 import var2
-from Wpyr import Wpyr
-from zconv2 import zconv2
+from .binomialFilter import binomialFilter
+from .blurDn import blurDn
+from .blur import blur
+from .cconv2 import cconv2
+from .clip import clip
+from .comparePyr import comparePyr
+from .compareRecon import compareRecon
+from .corrDn import corrDn
+from .entropy2 import entropy2
+from .factorial import factorial
+from .Gpyr import Gpyr
+from .histoMatch import histoMatch
+from .histo import histo
+from .idx2LB import idx2LB
+from .imGradient import imGradient
+from .imStats import imStats
+from .kurt2 import kurt2
+from .LB2idx import LB2idx
+from .Lpyr import Lpyr
+from .maxPyrHt import maxPyrHt
+from .mkAngle import mkAngle
+from .mkAngularSine import mkAngularSine
+from .mkDisc import mkDisc
+from .mkFract import mkFract
+from .mkGaussian import mkGaussian
+from .mkImpulse import mkImpulse
+from .mkRamp import mkRamp
+from .mkR import mkR
+from .mkSine import mkSine
+from .mkSquare import mkSquare
+from .mkZonePlate import mkZonePlate
+from .modulateFlip import modulateFlip
+from .namedFilter import namedFilter
+from .nextSz import nextSz
+from .pointOp import pointOp
+from .pyramid import pyramid
+from .range2 import range2
+from .rconv2 import rconv2
+from .rcosFn import rcosFn
+from .round import round
+from .roundVal import roundVal
+from .SCFpyr import SCFpyr
+from .SFpyr import SFpyr
+from .shift import shift
+from .showIm import showIm
+from .skew2 import skew2
+from .sp0Filters import sp0Filters
+from .sp1Filters import sp1Filters
+from .sp3Filters import sp3Filters
+from .sp5Filters import sp5Filters
+from .Spyr import Spyr
+# from SpyrMulti import SpyrMulti
+from .steer2HarmMtx import steer2HarmMtx
+from .steer import steer
+from .strictly_decreasing import strictly_decreasing
+from .upBlur import upBlur
+from .upConv import upConv
+from .var2 import var2
+from .Wpyr import Wpyr
+from .zconv2 import zconv2
 import ctypes
 import os
-import JBhelpers
+from . import JBhelpers
 
 # load the C library
 lib = ctypes.cdll.LoadLibrary(os.path.dirname(os.path.realpath(__file__)) +
diff --git a/__pycache__/Gpyr.cpython-35.pyc b/__pycache__/Gpyr.cpython-35.pyc
new file mode 100644
index 0000000..ba11661
Binary files /dev/null and b/__pycache__/Gpyr.cpython-35.pyc differ
diff --git a/__pycache__/JBhelpers.cpython-35.pyc b/__pycache__/JBhelpers.cpython-35.pyc
new file mode 100644
index 0000000..a7d58be
Binary files /dev/null and b/__pycache__/JBhelpers.cpython-35.pyc differ
diff --git a/__pycache__/LB2idx.cpython-35.pyc b/__pycache__/LB2idx.cpython-35.pyc
new file mode 100644
index 0000000..00f6ad7
Binary files /dev/null and b/__pycache__/LB2idx.cpython-35.pyc differ
diff --git a/__pycache__/Lpyr.cpython-35.pyc b/__pycache__/Lpyr.cpython-35.pyc
new file mode 100644
index 0000000..a496e72
Binary files /dev/null and b/__pycache__/Lpyr.cpython-35.pyc differ
diff --git a/__pycache__/SCFpyr.cpython-35.pyc b/__pycache__/SCFpyr.cpython-35.pyc
new file mode 100644
index 0000000..8a8de32
Binary files /dev/null and b/__pycache__/SCFpyr.cpython-35.pyc differ
diff --git a/__pycache__/SFpyr.cpython-35.pyc b/__pycache__/SFpyr.cpython-35.pyc
new file mode 100644
index 0000000..e6a3d69
Binary files /dev/null and b/__pycache__/SFpyr.cpython-35.pyc differ
diff --git a/__pycache__/Spyr.cpython-35.pyc b/__pycache__/Spyr.cpython-35.pyc
new file mode 100644
index 0000000..61afe1c
Binary files /dev/null and b/__pycache__/Spyr.cpython-35.pyc differ
diff --git a/__pycache__/Wpyr.cpython-35.pyc b/__pycache__/Wpyr.cpython-35.pyc
new file mode 100644
index 0000000..678cb5a
Binary files /dev/null and b/__pycache__/Wpyr.cpython-35.pyc differ
diff --git a/__pycache__/__init__.cpython-35.pyc b/__pycache__/__init__.cpython-35.pyc
new file mode 100644
index 0000000..49f0f15
Binary files /dev/null and b/__pycache__/__init__.cpython-35.pyc differ
diff --git a/__pycache__/binomialFilter.cpython-35.pyc b/__pycache__/binomialFilter.cpython-35.pyc
new file mode 100644
index 0000000..4ead7eb
Binary files /dev/null and b/__pycache__/binomialFilter.cpython-35.pyc differ
diff --git a/__pycache__/blur.cpython-35.pyc b/__pycache__/blur.cpython-35.pyc
new file mode 100644
index 0000000..c2f4647
Binary files /dev/null and b/__pycache__/blur.cpython-35.pyc differ
diff --git a/__pycache__/blurDn.cpython-35.pyc b/__pycache__/blurDn.cpython-35.pyc
new file mode 100644
index 0000000..1c50a55
Binary files /dev/null and b/__pycache__/blurDn.cpython-35.pyc differ
diff --git a/__pycache__/cconv2.cpython-35.pyc b/__pycache__/cconv2.cpython-35.pyc
new file mode 100644
index 0000000..0e8cca6
Binary files /dev/null and b/__pycache__/cconv2.cpython-35.pyc differ
diff --git a/__pycache__/clip.cpython-35.pyc b/__pycache__/clip.cpython-35.pyc
new file mode 100644
index 0000000..eb774c5
Binary files /dev/null and b/__pycache__/clip.cpython-35.pyc differ
diff --git a/__pycache__/comparePyr.cpython-35.pyc b/__pycache__/comparePyr.cpython-35.pyc
new file mode 100644
index 0000000..21e4276
Binary files /dev/null and b/__pycache__/comparePyr.cpython-35.pyc differ
diff --git a/__pycache__/compareRecon.cpython-35.pyc b/__pycache__/compareRecon.cpython-35.pyc
new file mode 100644
index 0000000..01ecb88
Binary files /dev/null and b/__pycache__/compareRecon.cpython-35.pyc differ
diff --git a/__pycache__/corrDn.cpython-35.pyc b/__pycache__/corrDn.cpython-35.pyc
new file mode 100644
index 0000000..22cbc54
Binary files /dev/null and b/__pycache__/corrDn.cpython-35.pyc differ
diff --git a/__pycache__/entropy2.cpython-35.pyc b/__pycache__/entropy2.cpython-35.pyc
new file mode 100644
index 0000000..adc9df4
Binary files /dev/null and b/__pycache__/entropy2.cpython-35.pyc differ
diff --git a/__pycache__/factorial.cpython-35.pyc b/__pycache__/factorial.cpython-35.pyc
new file mode 100644
index 0000000..6ed73a6
Binary files /dev/null and b/__pycache__/factorial.cpython-35.pyc differ
diff --git a/__pycache__/histo.cpython-35.pyc b/__pycache__/histo.cpython-35.pyc
new file mode 100644
index 0000000..5eea9dc
Binary files /dev/null and b/__pycache__/histo.cpython-35.pyc differ
diff --git a/__pycache__/histoMatch.cpython-35.pyc b/__pycache__/histoMatch.cpython-35.pyc
new file mode 100644
index 0000000..77cd4b5
Binary files /dev/null and b/__pycache__/histoMatch.cpython-35.pyc differ
diff --git a/__pycache__/idx2LB.cpython-35.pyc b/__pycache__/idx2LB.cpython-35.pyc
new file mode 100644
index 0000000..eadea8b
Binary files /dev/null and b/__pycache__/idx2LB.cpython-35.pyc differ
diff --git a/__pycache__/imGradient.cpython-35.pyc b/__pycache__/imGradient.cpython-35.pyc
new file mode 100644
index 0000000..31a7fe2
Binary files /dev/null and b/__pycache__/imGradient.cpython-35.pyc differ
diff --git a/__pycache__/imStats.cpython-35.pyc b/__pycache__/imStats.cpython-35.pyc
new file mode 100644
index 0000000..6f047ca
Binary files /dev/null and b/__pycache__/imStats.cpython-35.pyc differ
diff --git a/__pycache__/kurt2.cpython-35.pyc b/__pycache__/kurt2.cpython-35.pyc
new file mode 100644
index 0000000..701f1cc
Binary files /dev/null and b/__pycache__/kurt2.cpython-35.pyc differ
diff --git a/__pycache__/maxPyrHt.cpython-35.pyc b/__pycache__/maxPyrHt.cpython-35.pyc
new file mode 100644
index 0000000..eea28af
Binary files /dev/null and b/__pycache__/maxPyrHt.cpython-35.pyc differ
diff --git a/__pycache__/mkAngle.cpython-35.pyc b/__pycache__/mkAngle.cpython-35.pyc
new file mode 100644
index 0000000..338623d
Binary files /dev/null and b/__pycache__/mkAngle.cpython-35.pyc differ
diff --git a/__pycache__/mkAngularSine.cpython-35.pyc b/__pycache__/mkAngularSine.cpython-35.pyc
new file mode 100644
index 0000000..4b27711
Binary files /dev/null and b/__pycache__/mkAngularSine.cpython-35.pyc differ
diff --git a/__pycache__/mkDisc.cpython-35.pyc b/__pycache__/mkDisc.cpython-35.pyc
new file mode 100644
index 0000000..e5bf9dc
Binary files /dev/null and b/__pycache__/mkDisc.cpython-35.pyc differ
diff --git a/__pycache__/mkFract.cpython-35.pyc b/__pycache__/mkFract.cpython-35.pyc
new file mode 100644
index 0000000..81510f4
Binary files /dev/null and b/__pycache__/mkFract.cpython-35.pyc differ
diff --git a/__pycache__/mkGaussian.cpython-35.pyc b/__pycache__/mkGaussian.cpython-35.pyc
new file mode 100644
index 0000000..1ed2def
Binary files /dev/null and b/__pycache__/mkGaussian.cpython-35.pyc differ
diff --git a/__pycache__/mkImpulse.cpython-35.pyc b/__pycache__/mkImpulse.cpython-35.pyc
new file mode 100644
index 0000000..4055bc0
Binary files /dev/null and b/__pycache__/mkImpulse.cpython-35.pyc differ
diff --git a/__pycache__/mkR.cpython-35.pyc b/__pycache__/mkR.cpython-35.pyc
new file mode 100644
index 0000000..a0d2c62
Binary files /dev/null and b/__pycache__/mkR.cpython-35.pyc differ
diff --git a/__pycache__/mkRamp.cpython-35.pyc b/__pycache__/mkRamp.cpython-35.pyc
new file mode 100644
index 0000000..d1c4239
Binary files /dev/null and b/__pycache__/mkRamp.cpython-35.pyc differ
diff --git a/__pycache__/mkSine.cpython-35.pyc b/__pycache__/mkSine.cpython-35.pyc
new file mode 100644
index 0000000..62e99c8
Binary files /dev/null and b/__pycache__/mkSine.cpython-35.pyc differ
diff --git a/__pycache__/mkSquare.cpython-35.pyc b/__pycache__/mkSquare.cpython-35.pyc
new file mode 100644
index 0000000..3f2c473
Binary files /dev/null and b/__pycache__/mkSquare.cpython-35.pyc differ
diff --git a/__pycache__/mkZonePlate.cpython-35.pyc b/__pycache__/mkZonePlate.cpython-35.pyc
new file mode 100644
index 0000000..c767fe2
Binary files /dev/null and b/__pycache__/mkZonePlate.cpython-35.pyc differ
diff --git a/__pycache__/modulateFlip.cpython-35.pyc b/__pycache__/modulateFlip.cpython-35.pyc
new file mode 100644
index 0000000..21d073e
Binary files /dev/null and b/__pycache__/modulateFlip.cpython-35.pyc differ
diff --git a/__pycache__/namedFilter.cpython-35.pyc b/__pycache__/namedFilter.cpython-35.pyc
new file mode 100644
index 0000000..2f12065
Binary files /dev/null and b/__pycache__/namedFilter.cpython-35.pyc differ
diff --git a/__pycache__/nextSz.cpython-35.pyc b/__pycache__/nextSz.cpython-35.pyc
new file mode 100644
index 0000000..758a13b
Binary files /dev/null and b/__pycache__/nextSz.cpython-35.pyc differ
diff --git a/__pycache__/pointOp.cpython-35.pyc b/__pycache__/pointOp.cpython-35.pyc
new file mode 100644
index 0000000..39abe0e
Binary files /dev/null and b/__pycache__/pointOp.cpython-35.pyc differ
diff --git a/__pycache__/pyramid.cpython-35.pyc b/__pycache__/pyramid.cpython-35.pyc
new file mode 100644
index 0000000..95679fc
Binary files /dev/null and b/__pycache__/pyramid.cpython-35.pyc differ
diff --git a/__pycache__/range2.cpython-35.pyc b/__pycache__/range2.cpython-35.pyc
new file mode 100644
index 0000000..a089c28
Binary files /dev/null and b/__pycache__/range2.cpython-35.pyc differ
diff --git a/__pycache__/rconv2.cpython-35.pyc b/__pycache__/rconv2.cpython-35.pyc
new file mode 100644
index 0000000..b50d18c
Binary files /dev/null and b/__pycache__/rconv2.cpython-35.pyc differ
diff --git a/__pycache__/rcosFn.cpython-35.pyc b/__pycache__/rcosFn.cpython-35.pyc
new file mode 100644
index 0000000..6c2b1cc
Binary files /dev/null and b/__pycache__/rcosFn.cpython-35.pyc differ
diff --git a/__pycache__/round.cpython-35.pyc b/__pycache__/round.cpython-35.pyc
new file mode 100644
index 0000000..3ab2915
Binary files /dev/null and b/__pycache__/round.cpython-35.pyc differ
diff --git a/__pycache__/roundVal.cpython-35.pyc b/__pycache__/roundVal.cpython-35.pyc
new file mode 100644
index 0000000..4f6f2b9
Binary files /dev/null and b/__pycache__/roundVal.cpython-35.pyc differ
diff --git a/__pycache__/shift.cpython-35.pyc b/__pycache__/shift.cpython-35.pyc
new file mode 100644
index 0000000..e89baba
Binary files /dev/null and b/__pycache__/shift.cpython-35.pyc differ
diff --git a/__pycache__/showIm.cpython-35.pyc b/__pycache__/showIm.cpython-35.pyc
new file mode 100644
index 0000000..dd4048e
Binary files /dev/null and b/__pycache__/showIm.cpython-35.pyc differ
diff --git a/__pycache__/skew2.cpython-35.pyc b/__pycache__/skew2.cpython-35.pyc
new file mode 100644
index 0000000..c2eccb1
Binary files /dev/null and b/__pycache__/skew2.cpython-35.pyc differ
diff --git a/__pycache__/sp0Filters.cpython-35.pyc b/__pycache__/sp0Filters.cpython-35.pyc
new file mode 100644
index 0000000..a41ce6b
Binary files /dev/null and b/__pycache__/sp0Filters.cpython-35.pyc differ
diff --git a/__pycache__/sp1Filters.cpython-35.pyc b/__pycache__/sp1Filters.cpython-35.pyc
new file mode 100644
index 0000000..7e5a3fb
Binary files /dev/null and b/__pycache__/sp1Filters.cpython-35.pyc differ
diff --git a/__pycache__/sp3Filters.cpython-35.pyc b/__pycache__/sp3Filters.cpython-35.pyc
new file mode 100644
index 0000000..c1436cf
Binary files /dev/null and b/__pycache__/sp3Filters.cpython-35.pyc differ
diff --git a/__pycache__/sp5Filters.cpython-35.pyc b/__pycache__/sp5Filters.cpython-35.pyc
new file mode 100644
index 0000000..4ad8773
Binary files /dev/null and b/__pycache__/sp5Filters.cpython-35.pyc differ
diff --git a/__pycache__/steer.cpython-35.pyc b/__pycache__/steer.cpython-35.pyc
new file mode 100644
index 0000000..638894c
Binary files /dev/null and b/__pycache__/steer.cpython-35.pyc differ
diff --git a/__pycache__/steer2HarmMtx.cpython-35.pyc b/__pycache__/steer2HarmMtx.cpython-35.pyc
new file mode 100644
index 0000000..a1f31dc
Binary files /dev/null and b/__pycache__/steer2HarmMtx.cpython-35.pyc differ
diff --git a/__pycache__/strictly_decreasing.cpython-35.pyc b/__pycache__/strictly_decreasing.cpython-35.pyc
new file mode 100644
index 0000000..4b8a4c4
Binary files /dev/null and b/__pycache__/strictly_decreasing.cpython-35.pyc differ
diff --git a/__pycache__/upBlur.cpython-35.pyc b/__pycache__/upBlur.cpython-35.pyc
new file mode 100644
index 0000000..5b11b12
Binary files /dev/null and b/__pycache__/upBlur.cpython-35.pyc differ
diff --git a/__pycache__/upConv.cpython-35.pyc b/__pycache__/upConv.cpython-35.pyc
new file mode 100644
index 0000000..02959b1
Binary files /dev/null and b/__pycache__/upConv.cpython-35.pyc differ
diff --git a/__pycache__/var2.cpython-35.pyc b/__pycache__/var2.cpython-35.pyc
new file mode 100644
index 0000000..e86cccc
Binary files /dev/null and b/__pycache__/var2.cpython-35.pyc differ
diff --git a/__pycache__/zconv2.cpython-35.pyc b/__pycache__/zconv2.cpython-35.pyc
new file mode 100644
index 0000000..01a05d9
Binary files /dev/null and b/__pycache__/zconv2.cpython-35.pyc differ
diff --git a/binomialFilter.py b/binomialFilter.py
index dfe41d0..bc4e879 100644
--- a/binomialFilter.py
+++ b/binomialFilter.py
@@ -4,7 +4,7 @@
 def binomialFilter(size):
     ''' returns a vector of binomial coefficients of order (size-1) '''
     if size < 2:
-        print "Error: size argument must be larger than 1"
+        print("Error: size argument must be larger than 1")
         exit(1)
     
     kernel = numpy.array([[0.5], [0.5]])
diff --git a/blur.py b/blur.py
index fd45fb1..df8f77f 100644
--- a/blur.py
+++ b/blur.py
@@ -1,7 +1,7 @@
 import numpy
-from namedFilter import namedFilter
-from corrDn import corrDn
-from upConv import upConv
+from .namedFilter import namedFilter
+from .corrDn import corrDn
+from .upConv import upConv
 
 def blur(*args):
     ''' RES = blur(IM, LEVELS, FILT)
@@ -18,8 +18,8 @@ def blur(*args):
 
     # REQUIRED ARG:
     if len(args) == 0:
-        print "blur(IM, LEVELS, FILT)"
-        print "first argument is required"
+        print("blur(IM, LEVELS, FILT)")
+        print("first argument is required")
         exit(1)
     else:
         im = numpy.array(args[0])
@@ -31,7 +31,7 @@ def blur(*args):
         nlevs = 1
 
     if len(args) > 2:
-        if isinstance(args[2], basestring):
+        if isinstance(args[2], str):
             filt = namedFilter(args[2])
         else:
             filt = numpy.array(args[2])
@@ -49,7 +49,7 @@ def blur(*args):
         if len(im.shape) == 1 or im.shape[0] == 1 or im.shape[1] == 1: 
             # 1D image
             if len(filt) == 2 and (numpy.asarray(filt.shape) != 1).any():
-                print 'Error: can not apply 2D filter to 1D signal'
+                print('Error: can not apply 2D filter to 1D signal')
                 return
             
             imIn = corrDn(im, filt, 'reflect1', len(im))
diff --git a/blurDn.py b/blurDn.py
index cfd1751..aca322b 100644
--- a/blurDn.py
+++ b/blurDn.py
@@ -1,6 +1,6 @@
 import numpy
-from namedFilter import namedFilter
-from corrDn import corrDn
+from .namedFilter import namedFilter
+from .corrDn import corrDn
 
 def blurDn(*args):
     ''' RES = blurDn(IM, LEVELS, FILT)
@@ -15,7 +15,7 @@ def blurDn(*args):
         function res = blurDn(im, nlevs, filt)  '''
 
     if len(args) == 0:
-        print "Error: image input parameter required."
+        print("Error: image input parameter required.")
         return
 
     im = numpy.array(args[0])
@@ -28,7 +28,7 @@ def blurDn(*args):
 
     if len(args) > 2:
         filt = args[2]
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = namedFilter(filt)
     else:
         filt = namedFilter('binom5')
@@ -48,7 +48,7 @@ def blurDn(*args):
             # 1D image
             if len(filt.shape) > 1 and (filt.shape[1]!=1 and filt.shape[2]!=1):
                 # >1D filter
-                print 'Error: Cannot apply 2D filter to 1D signal'
+                print('Error: Cannot apply 2D filter to 1D signal')
                 return
             # orient filter and image correctly
             if im.shape[0] == 1:
diff --git a/cconv2.py b/cconv2.py
index 3150f47..2d32819 100644
--- a/cconv2.py
+++ b/cconv2.py
@@ -18,9 +18,9 @@ def cconv2(*args):
         Python port by Rob Young, 8/15  '''
     
     if len(args) < 2:
-        print 'Error: cconv2 requires two input matrices!'
-        print 'Usage: cconv2(matrix1, matrix2, center)'
-        print 'where center parameter is optional'
+        print('Error: cconv2 requires two input matrices!')
+        print('Usage: cconv2(matrix1, matrix2, center)')
+        print('where center parameter is optional')
         return
     else:
         a = numpy.array(args[0])
@@ -38,7 +38,7 @@ def cconv2(*args):
         large = b
         small = a
     else:
-        print 'Error: one matrix must be larger than the other in both dimensions!'
+        print('Error: one matrix must be larger than the other in both dimensions!')
         return
     
     ly = large.shape[0]
diff --git a/clip.py b/clip.py
index aa4905a..3e003ab 100644
--- a/clip.py
+++ b/clip.py
@@ -10,8 +10,8 @@ def clip(*args):
         ported to Python by Rob Young, 8/15  '''
     
     if len(args) == 0 or len(args) > 3:
-        print 'Usage: clip(im, minVal or Range, maxVal)'
-        print 'first input parameter is required'
+        print('Usage: clip(im, minVal or Range, maxVal)')
+        print('first input parameter is required')
         return
         
     im = numpy.array(args[0])
@@ -20,7 +20,7 @@ def clip(*args):
         minVal = 0;
         maxVal = 1;
     elif len(args) == 2:
-        if isinstance(args[1], (int, long, float)):
+        if isinstance(args[1], (int, float)):
             minVal = args[1]
             maxVal = args[1]+1
         else:
@@ -31,7 +31,7 @@ def clip(*args):
         maxVal = args[2]
         
     if maxVal < minVal:
-        print 'Error: maxVal cannot be less than minVal!'
+        print('Error: maxVal cannot be less than minVal!')
         return
 
     im[numpy.where(im < minVal)] = minVal
diff --git a/comparePyr.py b/comparePyr.py
index 0434078..4d619e8 100644
--- a/comparePyr.py
+++ b/comparePyr.py
@@ -19,7 +19,7 @@ def comparePyr(matPyr, pyPyr):
             pySz += sz[0] * sz[1]
 
     if(matSz != pySz):
-        print "size difference: %d != %d, returning 0" % (matSz, pySz)
+        print("size difference: %d != %d, returning 0" % (matSz, pySz))
         return 0
 
     # values are the same?
@@ -34,14 +34,14 @@ def comparePyr(matPyr, pyPyr):
         matTmp = numpy.reshape(matTmp, bandSz, order='F')
         matStart = matStart+matLen
         if (matTmp != pyPyr.pyr[idx]).any():
-            print "some pyramid elements not identical: checking..."
+            print("some pyramid elements not identical: checking...")
             for i in range(bandSz[0]):
                 for j in range(bandSz[1]):
                     if matTmp[i,j] != pyPyr.pyr[idx][i,j]:
                         if ( math.fabs(matTmp[i,j] - pyPyr.pyr[idx][i,j]) > 
                              prec ):
-                            print "failed level:%d element:%d %d value:%.15f %.15f" % (idx, i, j, matTmp[i,j], pyPyr.pyr[idx][i,j])
+                            print("failed level:%d element:%d %d value:%.15f %.15f" % (idx, i, j, matTmp[i,j], pyPyr.pyr[idx][i,j]))
                             return 0
-            print "same to at least %f" % prec
+            print("same to at least %f" % prec)
 
     return 1
diff --git a/compareRecon.py b/compareRecon.py
index 8961af4..69b62ad 100644
--- a/compareRecon.py
+++ b/compareRecon.py
@@ -8,15 +8,15 @@ def compareRecon(recon1, recon2):
     ## FIX: make precision a input parameter
     prec = -11   # desired precision
     if recon1.shape != recon2.shape:
-        print 'shape is different!'
-        print recon1.shape
-        print recon2.shape
+        print('shape is different!')
+        print(recon1.shape)
+        print(recon2.shape)
         return 0
 
     for i in range(recon1.shape[0]):
         for j in range(recon2.shape[1]):
             if numpy.absolute(recon1[i,j].real - recon2[i,j].real) > math.pow(10,-11):
-                print "real: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].real, recon2[i,j].real, numpy.absolute(recon1[i,j].real-recon2[i,j].real))
+                print("real: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].real, recon2[i,j].real, numpy.absolute(recon1[i,j].real-recon2[i,j].real)))
                 return 0
             ## FIX: need a better way to test
             # if we have many significant digits to the left of decimal we 
@@ -27,9 +27,9 @@ def compareRecon(recon1, recon2):
                     prec = prec + int(math.log(numpy.abs(recon1[i,j].imag), 10))
                     if prec > 0:
                         prec = -1
-                print prec
+                print(prec)
                 if numpy.absolute(recon1[i,j].imag - recon2[i,j].imag) > math.pow(10, prec):
-                    print "imag: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].imag, recon2[i,j].imag, numpy.absolute(recon1[i,j].imag-recon2[i,j].imag))
+                    print("imag: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].imag, recon2[i,j].imag, numpy.absolute(recon1[i,j].imag-recon2[i,j].imag)))
                     return 0
 
     return 1
diff --git a/convolve.h b/convolve.h
index 48d55f7..8f7b9d1 100755
--- a/convolve.h
+++ b/convolve.h
@@ -53,3 +53,6 @@ int internal_wrap_expand(image_type *image, image_type *filt, int x_fdim, int y_
 			 int x_start, int x_step, int x_stop, 
 			 int y_start, int y_step, int y_stop,
 			 image_type *result, int x_rdim, int y_rdim);
+void internal_pointop(double *im, double *res, int size, double *lut,
+		      int lutsize, double origin, double increment,
+		      int warnings);
diff --git a/corrDn.py b/corrDn.py
index 0b4ec80..456f763 100644
--- a/corrDn.py
+++ b/corrDn.py
@@ -10,7 +10,7 @@ def corrDn(image = None, filt = None, edges = 'reflect1', step = (1,1),
            start = (0,0), stop = None, result = None):
 
     if image is None or filt is None:
-        print 'Error: image and filter are required input parameters!'
+        print('Error: image and filter are required input parameters!')
         return
     else:
         image = image.copy()
@@ -23,8 +23,8 @@ def corrDn(image = None, filt = None, edges = 'reflect1', step = (1,1),
         stop = (image.shape[0], image.shape[1])
 
     if result is None:
-        rxsz = len(range(start[0], stop[0], step[0]))
-        rysz = len(range(start[1], stop[1], step[1]))
+        rxsz = len(list(range(start[0], stop[0], step[0])))
+        rysz = len(list(range(start[1], stop[1], step[1])))
         result = numpy.zeros((rxsz, rysz))
     else:
         result = numpy.array(result.copy())
diff --git a/factorial.py b/factorial.py
index a04a0c3..d3ea8a2 100644
--- a/factorial.py
+++ b/factorial.py
@@ -8,7 +8,7 @@ def factorial(*args):
         EPS, 11/02, Python port by Rob Young, 10/15  '''
 
     # if scalar input make it a single element array
-    if isinstance(args[0], (int, long, float)):
+    if isinstance(args[0], (int, float)):
         num = numpy.array([args[0]])
     else:
         num = numpy.array(args[0])
diff --git a/histo.py b/histo.py
index 56ab20b..18b1485 100644
--- a/histo.py
+++ b/histo.py
@@ -1,6 +1,6 @@
 import numpy
-from range2 import range2
-from round import round
+from .range2 import range2
+from .round import round
 
 def histo(*args):
     ''' [N,X] = histo(MTX, nbinsOrBinsize, binCenter);
@@ -25,8 +25,8 @@ def histo(*args):
         Eero Simoncelli, 3/97.  ported to Python by Rob Young, 8/15.  '''
 
     if len(args) == 0 or len(args) > 3:
-        print 'Usage: histo(mtx, nbins, binCtr)'
-        print 'first argument is required'
+        print('Usage: histo(mtx, nbins, binCtr)')
+        print('first argument is required')
         return
     else:
         mtx = args[0]
@@ -47,7 +47,7 @@ def histo(*args):
             tmpNbins = ( round(float(mx-binCtr) / float(binSize)) - 
                          round(float(mn-binCtr) / float(binSize)) )
             if tmpNbins != args[1]:
-                print 'Warning: Using %d bins instead of requested number (%d)' % (tmpNbins, args[1])
+                print('Warning: Using %d bins instead of requested number (%d)' % (tmpNbins, args[1]))
     else:
         binSize = float(mx-mn) / 101.0
 
@@ -60,7 +60,7 @@ def histo(*args):
     # numpy.histogram uses bin edges, not centers like Matlab's hist
     #bins = firstBin + binSize * numpy.array(range(tmpNbins+1))
     # compute bin edges
-    binsE = firstEdge + binSize * numpy.array(range(tmpNbins+1))
+    binsE = firstEdge + binSize * numpy.array(list(range(tmpNbins+1)))
     
     [N, X] = numpy.histogram(mtx, binsE)
 
diff --git a/imGradient.py b/imGradient.py
index 7143ff3..83166f8 100644
--- a/imGradient.py
+++ b/imGradient.py
@@ -1,5 +1,5 @@
 import numpy
-from corrDn import corrDn
+from .corrDn import corrDn
 
 def imGradient(*args):
     ''' [dx, dy] = imGradient(im, edges) 
@@ -18,8 +18,8 @@ def imGradient(*args):
         Python port by Rob Young, 10/15  '''
     
     if len(args) == 0 or len(args) > 2:
-        print 'Usage: imGradient(image, edges)'
-        print "'edges' argument is optional"
+        print('Usage: imGradient(image, edges)')
+        print("'edges' argument is optional")
     elif len(args) == 1:
         edges = "dont-compute"
     elif len(args) == 2:
diff --git a/imStats.py b/imStats.py
index 5e7d437..b8024bc 100644
--- a/imStats.py
+++ b/imStats.py
@@ -1,7 +1,7 @@
 import numpy
-from range2 import range2
-from var2 import var2
-from kurt2 import kurt2
+from .range2 import range2
+from .var2 import var2
+from .kurt2 import kurt2
 
 def imStats(*args):
     ''' Report image (matrix) statistics.
@@ -11,16 +11,16 @@ def imStats(*args):
         stdev of the difference, and also SNR (relative to IM1).  '''
 
     if len(args) == 0:
-        print 'Error: at least one input image is required'
+        print('Error: at least one input image is required')
         return
     elif len(args) == 1 and not numpy.isreal(args[0]).all():
-        print 'Error: input images must be real-valued matrices'
+        print('Error: input images must be real-valued matrices')
         return
     elif len(args) == 2 and ( not numpy.isreal(args[0]).all() or not numpy.isreal(args[1]).all()):
-        print 'Error: input images must be real-valued matrices'
+        print('Error: input images must be real-valued matrices')
         return
     elif len(args) > 2:
-        print 'Error: maximum of two input images allowed'
+        print('Error: maximum of two input images allowed')
         return
 
     if len(args) == 2:
@@ -32,17 +32,17 @@ def imStats(*args):
             snr = numpy.inf
         else:
             snr = 10 * numpy.log10(var2(args[0])/v)
-        print 'Difference statistics:'
-        print '  Range: [%d, %d]' % (mn, mx)
-        print '  Mean: %f,  Stdev (rmse): %f,  SNR (dB): %f' % (mean,
+        print('Difference statistics:')
+        print('  Range: [%d, %d]' % (mn, mx))
+        print('  Mean: %f,  Stdev (rmse): %f,  SNR (dB): %f' % (mean,
                                                                 numpy.sqrt(v),
-                                                                snr)
+                                                                snr))
     else:
         (mn, mx) = range2(args[0])
         mean = args[0].mean()
         var = var2(args[0])
         stdev = numpy.sqrt(var.real) + numpy.sqrt(var.imag)
         kurt = kurt2(args[0], mean, stdev**2)
-        print 'Image statistics:'
-        print '  Range: [%f, %f]' % (mn, mx)
-        print '  Mean: %f,  Stdev: %f,  Kurtosis: %f' % (mean, stdev, kurt)
+        print('Image statistics:')
+        print('  Range: [%f, %f]' % (mn, mx))
+        print('  Mean: %f,  Stdev: %f,  Kurtosis: %f' % (mean, stdev, kurt))
diff --git a/kurt2.py b/kurt2.py
index 4a40f84..8c77966 100644
--- a/kurt2.py
+++ b/kurt2.py
@@ -6,7 +6,7 @@ def kurt2(*args):
         MEAN (optional) and VAR (optional) make the computation faster.  '''
 
     if len(args) == 0:
-        print 'Error: input matrix is required'
+        print('Error: input matrix is required')
 
     if len(args) < 2:
         mn = args[0].mean()
diff --git a/makeChangelog.py b/makeChangelog.py
index 6ae0411..5316189 100755
--- a/makeChangelog.py
+++ b/makeChangelog.py
@@ -2,25 +2,25 @@
 # script makes a change log txt file from git logs
 # usage: makeChangelog.py > CHANGELOG.txt
 
-import commands
+import subprocess
 
-atags = commands.getoutput('git tag -n')
+atags = subprocess.getoutput('git tag -n')
 atags = atags.split('\n');
-tags = commands.getoutput('git tag -l')
+tags = subprocess.getoutput('git tag -l')
 tags = tags.split('\n');
 
 #for t in range(0, len(tags)-1):
 for t in range(len(tags)-1, -1, -1):
     if t == 0:
         #print "*** " + atags[t+1]
-        print "*** " + atags[t]
+        print("*** " + atags[t])
     else:
         #print '\n\n*** ' + atags[t+1]
-        print '\n\n*** ' + atags[t]
+        print('\n\n*** ' + atags[t])
     #commandStr = 'git log %s..%s --pretty=%s' % (tags[t], tags[t+1], '%s')
     commandStr = 'git log %s..%s --pretty=%s' % (tags[t-1], tags[t], '%s')
-    changes = commands.getoutput(commandStr)
+    changes = subprocess.getoutput(commandStr)
     changes = changes.split('\n')
     #changes = changes[::-1]
     for line in changes:
-        print '  + ' + line
+        print('  + ' + line)
diff --git a/mkAngle.py b/mkAngle.py
index bdbe675..29f9164 100644
--- a/mkAngle.py
+++ b/mkAngle.py
@@ -9,10 +9,10 @@ def mkAngle(*args):
     if len(args) > 0:
         sz = args[0]
         if not isinstance(sz, tuple):
-            sz = (sz, sz)
+            sz = tuple(sz)
     else:
-        print "Error: first input parameter 'size' is required!"
-        print "makeAngle(size, phase, origin)"
+        print("Error: first input parameter 'size' is required!")
+        print("makeAngle(size, phase, origin)")
         return
 
     # ------------------------------------------------------------
@@ -30,8 +30,8 @@ def mkAngle(*args):
 
     #------------------------------------------------------------------
 
-    (xramp, yramp) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-origin[1], 
-                                 (numpy.array(range(1,sz[0]+1)))-origin[0])
+    (xramp, yramp) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-origin[1], 
+                                 (numpy.array(list(range(1,sz[0]+1))))-origin[0])
     xramp = numpy.array(xramp)
     yramp = numpy.array(yramp)
 
diff --git a/mkAngularSine.py b/mkAngularSine.py
index b5b35cc..2e78bd8 100644
--- a/mkAngularSine.py
+++ b/mkAngularSine.py
@@ -1,5 +1,5 @@
 import numpy
-from mkAngle import mkAngle
+from .mkAngle import mkAngle
 
 def mkAngularSine(*args):
     ''' IM = mkAngularSine(SIZE, HARMONIC, AMPL, PHASE, ORIGIN)
@@ -13,15 +13,15 @@ def mkAngularSine(*args):
         Eero Simoncelli, 2/97.  Python port by Rob Young, 7/15.  '''
 
     if len(args) == 0:
-        print "mkAngularSine(SIZE, HARMONIC, AMPL, PHASE, ORIGIN)"
-        print "first argument is required"
+        print("mkAngularSine(SIZE, HARMONIC, AMPL, PHASE, ORIGIN)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
diff --git a/mkDisc.py b/mkDisc.py
index 8529ed1..bdadf67 100644
--- a/mkDisc.py
+++ b/mkDisc.py
@@ -1,8 +1,8 @@
 import numpy
 import sys
-from mkR import mkR
-from rcosFn import rcosFn
-from pointOp import pointOp
+from .mkR import mkR
+from .rcosFn import rcosFn
+from .pointOp import pointOp
 
 def mkDisc(*args):
     ''' IM = mkDisc(SIZE, RADIUS, ORIGIN, TWIDTH, VALS)
@@ -18,15 +18,15 @@ def mkDisc(*args):
         Eero Simoncelli, 6/96. Python port by Rob Young, 7/15.  '''
 
     if len(args) == 0:
-        print "mkDisc(SIZE, RADIUS, ORIGIN, TWIDTH, VALS)"
-        print "first argument is required"
+        print("mkDisc(SIZE, RADIUS, ORIGIN, TWIDTH, VALS)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
diff --git a/mkFract.py b/mkFract.py
index 0b28637..c3e163d 100644
--- a/mkFract.py
+++ b/mkFract.py
@@ -1,6 +1,6 @@
 import numpy
-from mkR import mkR
-from var2 import var2
+from .mkR import mkR
+from .var2 import var2
 
 def mkFract(*args):
     ''' Make a matrix of dimensions SIZE (a [Y X] 2-vector, or a scalar)
@@ -13,9 +13,9 @@ def mkFract(*args):
               Make this more efficient!   '''
 
     if len(args) == 0:
-        print 'Error: input parameter dims required'
+        print('Error: input parameter dims required')
     else:
-        if isinstance(args[0], (int, long)) or len(args[0]) == 1:
+        if isinstance(args[0], int) or len(args[0]) == 1:
             dims = (args[0], args[0])
         elif args[0] == 1:
             dims = (args[1], args[1])
@@ -42,7 +42,7 @@ def mkFract(*args):
     fres = numpy.fft.ifft2(fres)
 
     if abs(fres.imag).max() > 1e-10:
-        print 'Symmetry error in creating fractal'
+        print('Symmetry error in creating fractal')
     else:
         res = numpy.real(fres)
         res = res / numpy.sqrt(var2(res))
diff --git a/mkGaussian.py b/mkGaussian.py
index 922fa6e..84e89f2 100644
--- a/mkGaussian.py
+++ b/mkGaussian.py
@@ -14,15 +14,15 @@ def mkGaussian(*args):
         Eero Simoncelli, 6/96. Python port by Rob Young, 7/15.  '''
 
     if len(args) == 0:
-        print "mkRamp(SIZE, COVARIANCE, MEAN, AMPLITUDE)"
-        print "first argument is required"
+        print("mkRamp(SIZE, COVARIANCE, MEAN, AMPLITUDE)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -46,14 +46,14 @@ def mkGaussian(*args):
 
     #---------------------------------------------------------------
         
-    (xramp, yramp) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-mn[1], 
-                                    numpy.array(range(1,sz[0]+1))-mn[0])
+    (xramp, yramp) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-mn[1], 
+                                    numpy.array(list(range(1,sz[0]+1)))-mn[0])
 
-    if isinstance(cov, (int, long, float)):
+    if isinstance(cov, (int, float)):
         if 'norm' == ampl:
             ampl = 1.0 / (2.0 * numpy.pi * cov)
         e = ( (xramp**2) + (yramp**2) ) / ( -2.0 * cov )
-    elif len(cov) == 2 and isinstance(cov[0], (int, long, float)):
+    elif len(cov) == 2 and isinstance(cov[0], (int, float)):
         if 'norm' == ampl:
             if cov[0]*cov[1] < 0:
                 ampl = 1.0 / (2.0 * numpy.pi * 
diff --git a/mkImpulse.py b/mkImpulse.py
index d888fb6..9d4ac2e 100644
--- a/mkImpulse.py
+++ b/mkImpulse.py
@@ -4,8 +4,8 @@ def mkImpulse(*args):
     ''' create an image that is all zeros except for an impulse '''
 
     if(len(args) == 0):
-        print "mkImpulse(size, origin, amplitude)"
-        print "first input parameter is required"
+        print("mkImpulse(size, origin, amplitude)")
+        print("first input parameter is required")
         return
     
     if(isinstance(args[0], int)):
@@ -13,7 +13,7 @@ def mkImpulse(*args):
     elif(isinstance(args[0], tuple)):
         sz = args[0]
     else:
-        print "size parameter must be either an integer or a tuple"
+        print("size parameter must be either an integer or a tuple")
         return
 
     if(len(args) > 1):
diff --git a/mkR.py b/mkR.py
index 1f68802..4c5c372 100644
--- a/mkR.py
+++ b/mkR.py
@@ -8,12 +8,12 @@ def mkR(*args):
         Eero Simoncelli, 6/96.  Ported to Python by Rob Young, 5/14.  '''
     
     if len(args) == 0:
-        print 'Error: first input parameter is required!'
+        print('Error: first input parameter is required!')
         return
     else:
         sz = args[0]
 
-    if isinstance(sz, (int, long)) or len(sz) == 1:
+    if isinstance(sz, int) or len(sz) == 1:
         sz = (sz, sz)
 
     # -----------------------------------------------------------------
@@ -31,8 +31,8 @@ def mkR(*args):
 
     # -----------------------------------------------------------------
 
-    (xramp2, yramp2) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-origin[1], 
-                                   numpy.array(range(1,sz[0]+1))-origin[0])
+    (xramp2, yramp2) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-origin[1], 
+                                   numpy.array(list(range(1,sz[0]+1)))-origin[0])
 
     res = (xramp2**2 + yramp2**2)**(expt/2.0)
     
diff --git a/mkRamp.py b/mkRamp.py
index 12ebe38..746851b 100644
--- a/mkRamp.py
+++ b/mkRamp.py
@@ -11,15 +11,15 @@ def mkRamp(*args):
         optional '''
     
     if len(args) == 0:
-        print "mkRamp(SIZE, DIRECTION, SLOPE, INTERCEPT, ORIGIN)"
-        print "first argument is required"
+        print("mkRamp(SIZE, DIRECTION, SLOPE, INTERCEPT, ORIGIN)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -49,8 +49,8 @@ def mkRamp(*args):
     xinc = slope * math.cos(direction)
     yinc = slope * math.sin(direction)
 
-    [xramp, yramp] = numpy.meshgrid( xinc * (numpy.array(range(sz[1]))-origin[1]),
-                                  yinc * (numpy.array(range(sz[0]))-origin[0]) )
+    [xramp, yramp] = numpy.meshgrid( xinc * (numpy.array(list(range(sz[1])))-origin[1]),
+                                  yinc * (numpy.array(list(range(sz[0])))-origin[0]) )
 
     res = intercept + xramp + yramp
 
diff --git a/mkSine.py b/mkSine.py
index b6a5a80..d43ad20 100644
--- a/mkSine.py
+++ b/mkSine.py
@@ -1,6 +1,6 @@
 import numpy
 import math
-from mkRamp import mkRamp
+from .mkRamp import mkRamp
 
 def mkSine(*args):
     ''' IM = mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)
@@ -20,20 +20,20 @@ def mkSine(*args):
     # REQUIRED args:
 
     if len(args) < 2:
-        print "mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)"
-        print "       or"
-        print "mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN)"
-        print "first two arguments are required"
+        print("mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)")
+        print("       or")
+        print("mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN)")
+        print("first two arguments are required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
-    if isinstance(args[1], (int, float, long)):
+    if isinstance(args[1], (int, float)):
         frequency = (2.0 * numpy.pi) / args[1]
         # OPTIONAL args:
         if len(args) > 2:
diff --git a/mkSquare.py b/mkSquare.py
index 7c65554..4d09bbf 100644
--- a/mkSquare.py
+++ b/mkSquare.py
@@ -1,8 +1,8 @@
 import numpy
 import math
-from mkRamp import mkRamp
-from rcosFn import rcosFn
-from pointOp import pointOp
+from .mkRamp import mkRamp
+from .rcosFn import rcosFn
+from .pointOp import pointOp
 
 def mkSquare(*args):
     ''' IM = mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)
@@ -26,20 +26,20 @@ def mkSquare(*args):
     #REQUIRED ARGS:
 
     if len(args) < 2:
-        print "mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)"
-        print "       or"
-        print "mkSquare(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN, TWIDTH)"
-        print "first two arguments are required"
+        print("mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)")
+        print("       or")
+        print("mkSquare(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN, TWIDTH)")
+        print("first two arguments are required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
     
-    if isinstance(args[1], (int, float, long)):
+    if isinstance(args[1], (int, float)):
         frequency = (2.0 * numpy.pi) / args[1]
         # OPTIONAL args:
         if len(args) > 2:
diff --git a/mkZonePlate.py b/mkZonePlate.py
index 6040575..9c6663d 100644
--- a/mkZonePlate.py
+++ b/mkZonePlate.py
@@ -1,5 +1,5 @@
 import numpy
-from mkR import mkR
+from .mkR import mkR
 
 def mkZonePlate(*args):
     ''' IM = mkZonePlate(SIZE, AMPL, PHASE)
@@ -14,15 +14,15 @@ def mkZonePlate(*args):
     # REQUIRED ARGS:
 
     if len(args) == 0:
-        print "mkZonePlate(SIZE, AMPL, PHASE)"
-        print "first argument is required"
+        print("mkZonePlate(SIZE, AMPL, PHASE)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
     
     #---------------------------------------------------------------------
diff --git a/modulateFlip.py b/modulateFlip.py
index e2e4f75..d17367c 100644
--- a/modulateFlip.py
+++ b/modulateFlip.py
@@ -8,7 +8,7 @@ def modulateFlip(*args):
         (e.g., see Simoncelli90).  '''
 
     if len(args) == 0:
-        print "Error: filter input parameter required."
+        print("Error: filter input parameter required.")
         return
 
     lfilt = args[0]
@@ -18,13 +18,13 @@ def modulateFlip(*args):
     elif lfilt.shape[0] == 1:
         lfilt = lfilt.reshape(lfilt.shape[1], 1)
     elif len(lfilt.shape) > 2 or lfilt.shape[1] != 1:
-        print 'Error: only 1D input supported.'
+        print('Error: only 1D input supported.')
         return
 
     sz = len(lfilt)
     sz2 = numpy.ceil(sz/2.0);
 
-    ind = numpy.array(range(sz-1,-1,-1))
+    ind = numpy.array(list(range(sz-1,-1,-1)))
 
     hfilt = lfilt[ind].T * (-1)**((ind+1)-sz2)
 
diff --git a/namedFilter.py b/namedFilter.py
index 0b6e125..8af2894 100644
--- a/namedFilter.py
+++ b/namedFilter.py
@@ -1,4 +1,4 @@
-from binomialFilter import binomialFilter
+from .binomialFilter import binomialFilter
 import math
 import numpy
 
@@ -54,7 +54,7 @@ def namedFilter(name):
     elif name is "gauss3":  # for backward-compatibility
         kernel = math.sqrt(2) * numpy.array([[0.25], [0.5], [0.25]])
     else:
-        print "Error: Bad filter name: %s" % (name)
+        print("Error: Bad filter name: %s" % (name))
         exit(1)
 
     return numpy.array(kernel)
diff --git a/pointOp.py b/pointOp.py
index 8d9d4b9..c0025b7 100644
--- a/pointOp.py
+++ b/pointOp.py
@@ -9,12 +9,12 @@
 def pointOp(image, lut, origin, increment, warnings):
     result = numpy.zeros((image.shape[0], image.shape[1]))
 
-    lib.internal_pointop(image.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), 
+    lib.internal_pointop(image.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                          result.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
-                         image.shape[0] * image.shape[1], 
+                         image.shape[0] * image.shape[1],
                          lut.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
-                         lut.shape[0], 
-                         ctypes.c_double(origin), 
+                         lut.shape[0],
+                         ctypes.c_double(origin),
                          ctypes.c_double(increment), warnings)
 
     return numpy.array(result)
diff --git a/pyPyrTools.py b/pyPyrTools.py
index f1e5b71..ebb4fed 100644
--- a/pyPyrTools.py
+++ b/pyPyrTools.py
@@ -1,13 +1,13 @@
 import numpy
 #import pyPyrUtils as ppu
-import pyPyrUtils
+from . import pyPyrUtils
 #import pyPyrCcode
 import math
 import matplotlib.cm
 import os
 import scipy.misc
 import cmath
-import JBhelpers
+from . import JBhelpers
 import pylab
 import copy
 
@@ -20,7 +20,7 @@ class pyramid:  # pyramid
 
     # constructor
     def __init__(self):
-        print "please specify type of pyramid to create (Gpry, Lpyr, etc.)"
+        print("please specify type of pyramid to create (Gpry, Lpyr, etc.)")
         return
 
     # methods
@@ -41,7 +41,7 @@ def __init__(self, *args):    # (image height, filter file, edges)
         if len(args) > 0:
             self.image = numpy.array(args[0])
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -57,10 +57,10 @@ def __init__(self, *args):    # (image height, filter file, edges)
             elif args[2] == 'sp5Filters':
                 filters = pyPyrUtils.sp5Filters()
             elif os.path.isfile(args[2]):
-                print "Filter files not supported yet"
+                print("Filter files not supported yet")
                 return
             else:
-                print "filter parameters value %s not supported" % (args[2])
+                print("filter parameters value %s not supported" % (args[2]))
                 return
         else:
             filters = pyPyrUtils.sp1Filters()
@@ -77,8 +77,8 @@ def __init__(self, *args):    # (image height, filter file, edges)
             if args[1] == 'auto':
                 ht = max_ht
             elif args[1] > max_ht:
-                print "Error: cannot build pyramid higher than %d levels." % (
-                    max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (
+                    max_ht))
                 return
             else:
                 ht = args[1]
@@ -134,20 +134,20 @@ def __init__(self, *args):    # (image height, filter file, edges)
     # methods
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        if isinstance(args[1], (int, long)):
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        if isinstance(args[1], int):
             self.pyr[args[0]][0][args[1]] = args[2]
         elif isinstance(args[1], tuple):
             self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
         else:
-            print 'Error: location parameter must be int or tuple!'
+            print('Error: location parameter must be int or tuple!')
             return
 
     def spyrLev(self, lev):
         if lev < 0 or lev > self.spyrHt()-1:
-            print 'Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1)
+            print('Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1))
             return
         
         levArray = []
@@ -159,10 +159,10 @@ def spyrLev(self, lev):
 
     def spyrBand(self, lev, band):
         if lev < 0 or lev > self.spyrHt()-1:
-            print 'Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1)
+            print('Error: level parameter must be between 0 and %d!' % (self.spyrHt()-1))
             return
         if band < 0 or band > self.numBands()-1:
-            print 'Error: band parameter must be between 0 and %d!' % (self.numBands()-1)
+            print('Error: band parameter must be between 0 and %d!' % (self.numBands()-1))
 
         return self.band( ((lev*self.numBands())+band)+1 )
 
@@ -203,10 +203,10 @@ def reconPyr(self, *args):
             elif args[0] == 'sp5Filters':
                 filters = pyPyrUtils.sp5Filters()
             elif os.path.isfile(args[0]):
-                print "Filter files not supported yet"
+                print("Filter files not supported yet")
                 return
             else:
-                print "filter %s not supported" % (args[0])
+                print("filter %s not supported" % (args[0]))
                 return
         else:
             filters = pyPyrUtils.sp1Filters()
@@ -238,22 +238,22 @@ def reconPyr(self, *args):
         
         maxLev = 2 + self.spyrHt()
         if levs == 'all':
-            levs = numpy.array(range(maxLev))
+            levs = numpy.array(list(range(maxLev)))
         else:
             levs = numpy.array(levs)
             if (levs < 0).any() or (levs >= maxLev).any():
-                print "Error: level numbers must be in the range [0, %d]." % (maxLev-1)
+                print("Error: level numbers must be in the range [0, %d]." % (maxLev-1))
                 return
             else:
                 levs = numpy.array(levs)
                 if len(levs) > 1 and levs[0] < levs[1]:
                     levs = levs[::-1]  # we want smallest first
         if bands == 'all':
-            bands = numpy.array(range(self.numBands()))
+            bands = numpy.array(list(range(self.numBands())))
         else:
             bands = numpy.array(bands)
             if (bands < 0).any() or (bands > bfilts.shape[1]).any():
-                print "Error: band numbers must be in the range [0, %d]." % (self.numBands()-1)
+                print("Error: band numbers must be in the range [0, %d]." % (self.numBands()-1))
                 return
             else:
                 bands = numpy.array(bands)
@@ -383,10 +383,10 @@ def showPyr(self, prange = 'auto2', gap = 1, scale = 2, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.sqrt( numpy.var(band) )
             prange[nind-1,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(prange, basestring):
-            print "Error:Bad RANGE argument: %s'" % (prange)
+        elif isinstance(prange, str):
+            print("Error:Bad RANGE argument: %s'" % (prange))
         elif prange.shape[0] == 1 and prange.shape[1] == 2:
-            scales = numpy.power(scale, range(ht))
+            scales = numpy.power(scale, list(range(ht)))
             scales = numpy.outer( numpy.ones((nbands,1)), scales )
             scales = numpy.array([1, scales, numpy.power(scale, ht)])
             prange = numpy.outer(scales, prange)
@@ -405,11 +405,11 @@ def showPyr(self, prange = 'auto2', gap = 1, scale = 2, disp = 'qt'):
             ncols = int(numpy.ceil((nbands+1)/2))
             nrows = int(numpy.ceil(nbands/2))
 
-        a = numpy.array(range(1-nrows, 1))
+        a = numpy.array(list(range(1-nrows, 1)))
         b = numpy.zeros((1,ncols))[0]
         ab = numpy.concatenate((a,b))
         c = numpy.zeros((1,nrows))[0]
-        d = range(-1, -ncols-1, -1)
+        d = list(range(-1, -ncols-1, -1))
         cd = numpy.concatenate((c,d))
         relpos = numpy.vstack((ab,cd)).T
         
@@ -462,7 +462,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         if len(args) > 0:
             self.image = args[0]
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -471,7 +471,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         max_ht = numpy.floor( numpy.log2( min(self.image.shape) ) ) - 2
         if len(args) > 1:
             if(args[1] > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (max_ht))
             ht = args[1]
         else:
             ht = max_ht
@@ -479,7 +479,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             
         if len(args) > 2:
             if args[2] > 15 or args[2] < 0:
-                print "Warning: order must be an integer in the range [0,15]. Truncating."
+                print("Warning: order must be an integer in the range [0,15]. Truncating.")
                 order = min( max(args[2],0), 15 )
             else:
                 order = args[2]
@@ -490,7 +490,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
 
         if len(args) > 3:
             if args[3] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
+                print("Warning: twidth must be positive. Setting to 1.")
                 twidth = 1
             else:
                 twidth = args[3]
@@ -501,21 +501,21 @@ def __init__(self, *args):    # (image, height, order, twidth)
         # steering stuff:
 
         if nbands % 2 == 0:
-            harmonics = numpy.array(range(nbands/2)) * 2 + 1
+            harmonics = numpy.array(list(range(nbands/2))) * 2 + 1
         else:
-            harmonics = numpy.array(range((nbands-1)/2)) * 2
+            harmonics = numpy.array(list(range((nbands-1)/2))) * 2
 
         steermtx = pyPyrUtils.steer2HarmMtx(harmonics, 
-                                            numpy.pi*numpy.array(range(nbands))/nbands,
+                                            numpy.pi*numpy.array(list(range(nbands)))/nbands,
                                             'even')
         #------------------------------------------------------
         
         dims = numpy.array(self.image.shape)
         ctr = numpy.ceil((numpy.array(dims)+0.5)/2)
         
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
+        (xramp, yramp) = numpy.meshgrid((numpy.array(list(range(1,dims[1]+1)))-ctr[1])/
                                      (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
+                                     (numpy.array(list(range(1,dims[0]+1)))-ctr[0])/
                                      (dims[0]/2))
         angle = numpy.arctan2(yramp, xramp)
         log_rad = numpy.sqrt(xramp**2 + yramp**2)
@@ -556,7 +556,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             Xrcos -= numpy.log2(2)
 
             lutsize = 1024
-            Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
+            Xcosn = numpy.pi * numpy.array(list(range(-(2*lutsize+1), (lutsize+2)))) / lutsize
 
             order = nbands -1
             const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
@@ -639,7 +639,7 @@ def reconSFpyr(self, *args):
 
         if len(args) > 2:
             if args[2] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
+                print("Warning: twidth must be positive. Setting to 1.")
                 twidth = 1
             else:
                 twidth = args[2]
@@ -651,20 +651,20 @@ def reconSFpyr(self, *args):
         nbands = self.numBands()
         
         maxLev = 1 + self.spyrHt()
-        if isinstance(levs, basestring) and levs == 'all':
-            levs = numpy.array(range(maxLev+1))
-        elif isinstance(levs, basestring):
-            print "Error: %s not valid for levs parameter." % (levs)
-            print "levs must be either a 1D numpy array or the string 'all'."
+        if isinstance(levs, str) and levs == 'all':
+            levs = numpy.array(list(range(maxLev+1)))
+        elif isinstance(levs, str):
+            print("Error: %s not valid for levs parameter." % (levs))
+            print("levs must be either a 1D numpy array or the string 'all'.")
             return
         else:
             levs = numpy.array(levs)
 
-        if isinstance(bands, basestring) and bands == 'all':
-            bands = numpy.array(range(nbands))
-        elif isinstance(bands, basestring):
-            print "Error: %s not valid for bands parameter." % (bands)
-            print "bands must be either a 1D numpy array or the string 'all'."
+        if isinstance(bands, str) and bands == 'all':
+            bands = numpy.array(list(range(nbands)))
+        elif isinstance(bands, str):
+            print("Error: %s not valid for bands parameter." % (bands))
+            print("bands must be either a 1D numpy array or the string 'all'.")
             return
         else:
             bands = numpy.array(bands)
@@ -693,9 +693,9 @@ def reconSFpyr(self, *args):
         dims = numpy.array(self.pyrSize[0])
         ctr = numpy.ceil((dims+0.5)/2.0)
 
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
+        (xramp, yramp) = numpy.meshgrid((numpy.array(list(range(1,dims[1]+1)))-ctr[1])/
                                      (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
+                                     (numpy.array(list(range(1,dims[0]+1)))-ctr[0])/
                                      (dims[0]/2))
         angle = numpy.arctan2(yramp, xramp)
         log_rad = numpy.sqrt(xramp**2 + yramp**2)
@@ -709,7 +709,7 @@ def reconSFpyr(self, *args):
 
         # from reconSFpyrLevs
         lutsize = 1024
-        Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
+        Xcosn = numpy.pi * numpy.array(list(range(-(2*lutsize+1), (lutsize+2)))) / lutsize
         
         order = nbands -1
         const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
@@ -852,7 +852,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         if len(args) > 0:
             self.image = args[0]
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -861,7 +861,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         max_ht = numpy.floor( numpy.log2( min(self.image.shape) ) ) - 2
         if len(args) > 1:
             if(args[1] > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (max_ht))
             ht = args[1]
         else:
             ht = max_ht
@@ -869,7 +869,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             
         if len(args) > 2:
             if args[2] > 15 or args[2] < 0:
-                print "Warning: order must be an integer in the range [0,15]. Truncating."
+                print("Warning: order must be an integer in the range [0,15]. Truncating.")
                 order = min( max(args[2],0), 15 )
             else:
                 order = args[2]
@@ -880,7 +880,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
 
         if len(args) > 3:
             if args[3] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
+                print("Warning: twidth must be positive. Setting to 1.")
                 twidth = 1
             else:
                 twidth = args[3]
@@ -891,21 +891,21 @@ def __init__(self, *args):    # (image, height, order, twidth)
         # steering stuff:
 
         if nbands % 2 == 0:
-            harmonics = numpy.array(range(nbands/2)) * 2 + 1
+            harmonics = numpy.array(list(range(nbands/2))) * 2 + 1
         else:
-            harmonics = numpy.array(range((nbands-1)/2)) * 2
+            harmonics = numpy.array(list(range((nbands-1)/2))) * 2
 
         steermtx = pyPyrUtils.steer2HarmMtx(harmonics, 
-                                     numpy.pi*numpy.array(range(nbands))/nbands,
+                                     numpy.pi*numpy.array(list(range(nbands)))/nbands,
                                      'even')
         #------------------------------------------------------
         
         dims = numpy.array(self.image.shape)
         ctr = numpy.ceil((numpy.array(dims)+0.5)/2)
         
-        (xramp, yramp) = numpy.meshgrid((numpy.array(range(1,dims[1]+1))-ctr[1])/
+        (xramp, yramp) = numpy.meshgrid((numpy.array(list(range(1,dims[1]+1)))-ctr[1])/
                                      (dims[1]/2), 
-                                     (numpy.array(range(1,dims[0]+1))-ctr[0])/
+                                     (numpy.array(list(range(1,dims[0]+1)))-ctr[0])/
                                      (dims[0]/2))
         angle = numpy.arctan2(yramp, xramp)
         log_rad = numpy.sqrt(xramp**2 + yramp**2)
@@ -946,7 +946,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             Xrcos -= numpy.log2(2)
 
             lutsize = 1024
-            Xcosn = numpy.pi * numpy.array(range(-(2*lutsize+1), (lutsize+2))) / lutsize
+            Xcosn = numpy.pi * numpy.array(list(range(-(2*lutsize+1), (lutsize+2)))) / lutsize
 
             order = nbands -1
             const = (2**(2*order))*(scipy.misc.factorial(order, exact=True)**2)/float(nbands*scipy.misc.factorial(2*order, exact=True))
@@ -1010,7 +1010,7 @@ def reconPyr(self, *args):
 
         if len(args) > 2:
             if args[2] <= 0:
-                print "Warning: twidth must be positive. Setting to 1."
+                print("Warning: twidth must be positive. Setting to 1.")
                 twidth = 1
             else:
                 twidth = args[2]
@@ -1064,17 +1064,17 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
         if len(args) > 0:
             self.image = args[0]
         else:
-            print "pyr = Lpyr(image, height, filter1, filter2, edges)"
-            print "First argument (image) is required"
+            print("pyr = Lpyr(image, height, filter1, filter2, edges)")
+            print("First argument (image) is required")
             return
 
         if len(args) > 2:
             filt1 = args[2]
-            if isinstance(filt1, basestring):
+            if isinstance(filt1, str):
                 filt1 = pyPyrUtils.namedFilter(filt1)
             elif len(filt1.shape) != 1 and ( filt1.shape[0] != 1 and
                                              filt1.shape[1] != 1 ):
-                print "Error: filter1 should be a 1D filter (i.e., a vector)"
+                print("Error: filter1 should be a 1D filter (i.e., a vector)")
                 return
         else:
             filt1 = pyPyrUtils.namedFilter('binom5')
@@ -1085,11 +1085,11 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
 
         if len(args) > 3:
             filt2 = args[3]
-            if isinstance(filt2, basestring):
+            if isinstance(filt2, str):
                 filt2 = pyPyrUtils.namedFilter(filt2)
             elif len(filt2.shape) != 1 and ( filt2.shape[0] != 1 and
                                             filt2.shape[1] != 1 ):
-                print "Error: filter2 should be a 1D filter (i.e., a vector)"
+                print("Error: filter2 should be a 1D filter (i.e., a vector)")
                 return
         else:
             filt2 = filt1
@@ -1102,8 +1102,8 @@ def __init__(self, *args):    # (image, height, filter1, filter2, edges)
             else:
                 self.height = args[1]
                 if self.height > maxHeight:
-                    print ( "Error: cannot build pyramid higher than %d levels"
-                            % (maxHeight) )
+                    print(( "Error: cannot build pyramid higher than %d levels"
+                            % (maxHeight) ))
                     return
         else:
             self.height = maxHeight
@@ -1192,16 +1192,16 @@ def catBands(self, *args):
     # set a pyramid value
     def set_old(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, element, value)'
-        print 'band=%d  element=%d  value=%d' % (args[0],args[1],args[2])
-        print self.pyr[args[0]].shape
+            print('Error: three input parameters required:')
+            print('  set(band, element, value)')
+        print('band=%d  element=%d  value=%d' % (args[0],args[1],args[2]))
+        print(self.pyr[args[0]].shape)
         self.pyr[args[0]][args[1]] = args[2] 
 
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, element(tuple), value)'
+            print('Error: three input parameters required:')
+            print('  set(band, element(tuple), value)')
         self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
 
     def reconPyr(self, *args):
@@ -1223,14 +1223,14 @@ def reconPyr(self, *args):
         maxLev = self.height
 
         if levs == 'all':
-            levs = range(0,maxLev)
+            levs = list(range(0,maxLev))
         else:
             if (levs > maxLev-1).any():
-                print ( "Error: level numbers must be in the range [0, %d]." % 
-                        (maxLev-1) )
+                print(( "Error: level numbers must be in the range [0, %d]." % 
+                        (maxLev-1) ))
                 return
 
-        if isinstance(filt2, basestring):
+        if isinstance(filt2, str):
             filt2 = pyPyrUtils.namedFilter(filt2)
         else:
             if len(filt2.shape) == 1:
@@ -1354,11 +1354,11 @@ def showPyr(self, pRange = None, gap = 1, scale = None, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.std(band)
             pRange[nind,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(pRange, basestring):
-            print "Error: band range argument: %s" % (pRange)
+        elif isinstance(pRange, str):
+            print("Error: band range argument: %s" % (pRange))
             return
         elif pRange.shape[0] == 1 and pRange.shape[1] == 2:
-            scales = numpy.power( numpy.array( range(0,nind) ), scale)
+            scales = numpy.power( numpy.array( list(range(0,nind)) ), scale)
             pRange = numpy.outer( scales, pRange )
             band = self.pyrLow()
             pRange[nind,:] = ( pRange[nind,:] + numpy.mean(band) - 
@@ -1405,7 +1405,7 @@ def showPyr(self, pRange = None, gap = 1, scale = None, disp = 'qt'):
                 llpos[bnum,:] = ctr - numpy.floor(numpy.array(sz))/2.0 
             # make position list positive, and allocate appropriate image
             llpos = llpos - numpy.ones((nind,1))*numpy.min(llpos)
-            pind = range(self.height)
+            pind = list(range(self.height))
             for i in pind:
                 pind[i] = self.band(i).shape
             urpos = llpos + pind
@@ -1433,8 +1433,8 @@ class Gpyr(Lpyr):
     def __init__(self, *args):    # (image, height, filter, edges)
         self.pyrType = 'Gaussian'
         if len(args) < 1:
-            print "pyr = Gpyr(image, height, filter, edges)"
-            print "First argument (image) is required"
+            print("pyr = Gpyr(image, height, filter, edges)")
+            print("First argument (image) is required")
             return
         else:
             self.image = args[0]
@@ -1442,10 +1442,10 @@ def __init__(self, *args):    # (image, height, filter, edges)
         if len(args) > 2:
             filt = args[2]
             if not (filt.shape == 1).any():
-                print "Error: filt should be a 1D filter (i.e., a vector)"
+                print("Error: filt should be a 1D filter (i.e., a vector)")
                 return
         else:
-            print "no filter set, so filter is binom5"
+            print("no filter set, so filter is binom5")
             filt = pyPyrUtils.namedFilter('binom5')
             if self.image.shape[0] == 1:
                 filt = filt.reshape(1,5)
@@ -1460,8 +1460,8 @@ def __init__(self, *args):    # (image, height, filter, edges)
             else:
                 self.height = args[1]
                 if self.height > maxHeight:
-                    print ( "Error: cannot build pyramid higher than %d levels"
-                            % (maxHeight) )
+                    print(( "Error: cannot build pyramid higher than %d levels"
+                            % (maxHeight) ))
                     return
         else:
             self.height = maxHeight
@@ -1521,7 +1521,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
         if len(args) > 0:
             im = args[0]
         else:
-            print "First argument (image) is required."
+            print("First argument (image) is required.")
             return
 
         #------------------------------------------------
@@ -1531,11 +1531,11 @@ def __init__(self, *args):    # (image, height, order, twidth)
             filt = args[2]
         else:
             filt = "qmf9"
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = pyPyrUtils.namedFilter(filt)
 
         if len(filt.shape) != 1 and filt.shape[0] != 1 and filt.shape[1] != 1:
-            print "Error: filter should be 1D (i.e., a vector)";
+            print("Error: filter should be 1D (i.e., a vector)");
             return
         hfilt = pyPyrUtils.modulateFlip(filt)
         #hfilt = pyPyrUtils.modulateFlip(filt).T
@@ -1565,7 +1565,7 @@ def __init__(self, *args):    # (image, height, order, twidth)
             if ht == 'auto':
                 ht = max_ht
             elif(ht > max_ht):
-                print "Error: cannot build pyramid higher than %d levels." % (max_ht)
+                print("Error: cannot build pyramid higher than %d levels." % (max_ht))
         else:
             ht = max_ht
         ht = int(ht)
@@ -1678,28 +1678,28 @@ def reconPyr(self, *args):
         maxLev = int(self.wpyrHt() + 1)
 
         if levs == 'all':
-            levs = numpy.array(range(maxLev))
+            levs = numpy.array(list(range(maxLev)))
         else:
             tmpLevs = []
             for l in levs:
                 tmpLevs.append((maxLev-1)-l)
             levs = numpy.array(tmpLevs)
             if (levs > maxLev).any():
-                print "Error: level numbers must be in the range [0, %d]" % (maxLev)
-        allLevs = numpy.array(range(maxLev))
+                print("Error: level numbers must be in the range [0, %d]" % (maxLev))
+        allLevs = numpy.array(list(range(maxLev)))
 
         if bands == "all":
             if ( len(self.band(0)) == 1 or self.band(0).shape[0] == 1 or 
                  self.band(0).shape[1] == 1 ):
                 bands = numpy.array([0]);
             else:
-                bands = numpy.array(range(3))
+                bands = numpy.array(list(range(3)))
         else:
             bands = numpy.array(bands)
             if (bands < 0).any() or (bands > 2).any():
-                print "Error: band numbers must be in the range [0,2]."
+                print("Error: band numbers must be in the range [0,2].")
         
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = pyPyrUtils.namedFilter(filt)
 
         hfilt = pyPyrUtils.modulateFlip(filt).T
@@ -1858,25 +1858,25 @@ def reconPyr(self, *args):
 
     def set(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        if isinstance(args[1], (int, long)):
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        if isinstance(args[1], int):
             self.pyr[args[0]][0][args[1]] = args[2]
         elif isinstance(args[1], tuple):
             self.pyr[args[0]][args[1][0]][args[1][1]] = args[2] 
         else:
-            print 'Error: location parameter must be int or tuple!'
+            print('Error: location parameter must be int or tuple!')
             return
             
 
     def set1D(self, *args):
         if len(args) != 3:
-            print 'Error: three input parameters required:'
-            print '  set(band, location, value)'
-            print '  where band and value are integer and location is a tuple'
-        print '%d %d %d' % (args[0], args[1], args[2])
-        print self.pyr[args[0]][0][1]
+            print('Error: three input parameters required:')
+            print('  set(band, location, value)')
+            print('  where band and value are integer and location is a tuple')
+        print('%d %d %d' % (args[0], args[1], args[2]))
+        print(self.pyr[args[0]][0][1])
 
     def pyrLow(self):
         return numpy.array(self.band(len(self.pyrSize)-1))
@@ -1967,10 +1967,10 @@ def showPyr(self, prange = None, gap = 1, scale = None, disp = 'qt'):
             av = numpy.mean(band)
             stdev = numpy.sqrt( numpy.var(band) )
             prange[nind-1,:] = numpy.array([av-2*stdev, av+2*stdev])
-        elif isinstance(prange, basestring):
-            print "Error:Bad RANGE argument: %s'" % (prange)
+        elif isinstance(prange, str):
+            print("Error:Bad RANGE argument: %s'" % (prange))
         elif prange.shape[0] == 1 and prange.shape[1] == 2:
-            scales = numpy.power(scale, range(ht))
+            scales = numpy.power(scale, list(range(ht)))
             scales = numpy.outer( numpy.ones((nbands,1)), scales )
             scales = numpy.array([1, scales, numpy.power(scale, ht)])
             prange = numpy.outer(scales, prange)
diff --git a/pyPyrUtils.py b/pyPyrUtils.py
index 536d154..ee9c6ab 100644
--- a/pyPyrUtils.py
+++ b/pyPyrUtils.py
@@ -11,10 +11,10 @@
 import sys
 from PyQt4 import QtGui
 from PyQt4 import QtCore
-import JBhelpers
+from . import JBhelpers
 import PIL
-import ImageTk
-import Tkinter
+from PIL import ImageTk
+import tkinter
 import ctypes
 lib = ctypes.cdll.LoadLibrary('./wrapConv.so')
 
@@ -37,7 +37,7 @@ def maxPyrHt_old(imsz, filtsz):
         imsz = (imsz[0], 1)
         filtsz = (filtsz[0], 1)
     elif len(imsz) == 1 and not any(f == 1 for f in filtsz):
-            print "Error: cannot have a 1D 'image' and 2D filter"
+            print("Error: cannot have a 1D 'image' and 2D filter")
             exit(1)
     elif len(imsz) == 1:
         imsz = (imsz[0], 1)
@@ -109,7 +109,7 @@ def maxPyrHt(imsz, filtsz):
 # returns a vector of binomial coefficients of order (size-1)
 def binomialFilter(size):
     if size < 2:
-        print "Error: size argument must be larger than 1"
+        print("Error: size argument must be larger than 1")
         exit(1)
     
     kernel = numpy.array([[0.5], [0.5]])
@@ -171,7 +171,7 @@ def namedFilter(name):
     elif name is "gauss3":  # for backward-compatibility
         kernel = math.sqrt(2) * numpy.array([[0.25], [0.5], [0.25]])
     else:
-        print "Error: Bad filter name: %s" % (name)
+        print("Error: Bad filter name: %s" % (name))
         exit(1)
     return numpy.array(kernel)
 
@@ -181,15 +181,15 @@ def strictly_decreasing(L):
 def compareRecon(recon1, recon2):
     prec = -11
     if recon1.shape != recon2.shape:
-        print 'shape is different!'
-        print recon1.shape
-        print recon2.shape
+        print('shape is different!')
+        print(recon1.shape)
+        print(recon2.shape)
         return 0
 
     for i in range(recon1.shape[0]):
         for j in range(recon2.shape[1]):
             if numpy.absolute(recon1[i,j].real - recon2[i,j].real) > math.pow(10,-11):
-                print "real: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].real, recon2[i,j].real, numpy.absolute(recon1[i,j].real-recon2[i,j].real))
+                print("real: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].real, recon2[i,j].real, numpy.absolute(recon1[i,j].real-recon2[i,j].real)))
                 return 0
             ## FIX: need a better way to test
             # if we have many significant digits to the left of decimal we 
@@ -200,9 +200,9 @@ def compareRecon(recon1, recon2):
                     prec = prec + int(math.log(numpy.abs(recon1[i,j].imag), 10))
                     if prec > 0:
                         prec = -1
-                print prec
+                print(prec)
                 if numpy.absolute(recon1[i,j].imag - recon2[i,j].imag) > math.pow(10, prec):
-                    print "imag: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].imag, recon2[i,j].imag, numpy.absolute(recon1[i,j].imag-recon2[i,j].imag))
+                    print("imag: i=%d j=%d %.15f %.15f diff=%.15f" % (i, j, recon1[i,j].imag, recon2[i,j].imag, numpy.absolute(recon1[i,j].imag-recon2[i,j].imag)))
                     return 0
 
     return 1
@@ -221,7 +221,7 @@ def comparePyr(matPyr, pyPyr):
             pySz += sz[0] * sz[1]
 
     if(matSz != pySz):
-        print "size difference: %d != %d, returning 0" % (matSz, pySz)
+        print("size difference: %d != %d, returning 0" % (matSz, pySz))
         return 0
 
     # values are the same?
@@ -236,15 +236,15 @@ def comparePyr(matPyr, pyPyr):
         matTmp = numpy.reshape(matTmp, bandSz, order='F')
         matStart = matStart+matLen
         if (matTmp != pyPyr.pyr[idx]).any():
-            print "some pyramid elements not identical: checking..."
+            print("some pyramid elements not identical: checking...")
             for i in range(bandSz[0]):
                 for j in range(bandSz[1]):
                     if matTmp[i,j] != pyPyr.pyr[idx][i,j]:
                         if ( math.fabs(matTmp[i,j] - pyPyr.pyr[idx][i,j]) > 
                              prec ):
-                            print "failed level:%d element:%d %d value:%.15f %.15f" % (idx, i, j, matTmp[i,j], pyPyr.pyr[idx][i,j])
+                            print("failed level:%d element:%d %d value:%.15f %.15f" % (idx, i, j, matTmp[i,j], pyPyr.pyr[idx][i,j]))
                             return 0
-            print "same to at least %f" % prec
+            print("same to at least %f" % prec)
 
     return 1
 
@@ -260,15 +260,15 @@ def mkAngularSine(*args):
     # Eero Simoncelli, 2/97.  Python port by Rob Young, 7/15.
 
     if len(args) == 0:
-        print "mkAngularSine(SIZE, HARMONIC, AMPL, PHASE, ORIGIN)"
-        print "first argument is required"
+        print("mkAngularSine(SIZE, HARMONIC, AMPL, PHASE, ORIGIN)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -311,15 +311,15 @@ def mkGaussian(*args):
 # Eero Simoncelli, 6/96. Python port by Rob Young, 7/15.
 
     if len(args) == 0:
-        print "mkRamp(SIZE, COVARIANCE, MEAN, AMPLITUDE)"
-        print "first argument is required"
+        print("mkRamp(SIZE, COVARIANCE, MEAN, AMPLITUDE)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -343,14 +343,14 @@ def mkGaussian(*args):
 
     #---------------------------------------------------------------
         
-    (xramp, yramp) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-mn[1], 
-                                    numpy.array(range(1,sz[0]+1))-mn[0])
+    (xramp, yramp) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-mn[1], 
+                                    numpy.array(list(range(1,sz[0]+1)))-mn[0])
 
-    if isinstance(cov, (int, long, float)):
+    if isinstance(cov, (int, float)):
         if 'norm' == ampl:
             ampl = 1.0 / (2.0 * numpy.pi * cov)
         e = ( (xramp**2) + (yramp**2) ) / ( -2.0 * cov )
-    elif len(cov) == 2 and isinstance(cov[0], (int, long, float)):
+    elif len(cov) == 2 and isinstance(cov[0], (int, float)):
         if 'norm' == ampl:
             if cov[0]*cov[1] < 0:
                 ampl = 1.0 / (2.0 * numpy.pi * 
@@ -391,15 +391,15 @@ def mkDisc(*args):
  # Eero Simoncelli, 6/96. Python port by Rob Young, 7/15.
 
     if len(args) == 0:
-        print "mkDisc(SIZE, RADIUS, ORIGIN, TWIDTH, VALS)"
-        print "first argument is required"
+        print("mkDisc(SIZE, RADIUS, ORIGIN, TWIDTH, VALS)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -454,20 +454,20 @@ def mkSine(*args):
     # REQUIRED args:
 
     if len(args) < 2:
-        print "mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)"
-        print "       or"
-        print "mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN)"
-        print "first two arguments are required"
+        print("mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)")
+        print("       or")
+        print("mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN)")
+        print("first two arguments are required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
-    if isinstance(args[1], (int, float, long)):
+    if isinstance(args[1], (int, float)):
         frequency = (2.0 * numpy.pi) / args[1]
         # OPTIONAL args:
         if len(args) > 2:
@@ -526,15 +526,15 @@ def mkZonePlate(*args):
     # REQUIRED ARGS:
 
     if len(args) == 0:
-        print "mkZonePlate(SIZE, AMPL, PHASE)"
-        print "first argument is required"
+        print("mkZonePlate(SIZE, AMPL, PHASE)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
     
     #---------------------------------------------------------------------
@@ -576,20 +576,20 @@ def mkSquare(*args):
     # REQUIRED ARGS:
 
     if len(args) < 2:
-        print "mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)"
-        print "       or"
-        print "mkSquare(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN, TWIDTH)"
-        print "first two arguments are required"
+        print("mkSquare(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN, TWIDTH)")
+        print("       or")
+        print("mkSquare(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN, TWIDTH)")
+        print("first two arguments are required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
     
-    if isinstance(args[1], (int, float, long)):
+    if isinstance(args[1], (int, float)):
         frequency = (2.0 * numpy.pi) / args[1]
         # OPTIONAL args:
         if len(args) > 2:
@@ -659,15 +659,15 @@ def mkRamp(*args):
     # optional
     
     if len(args) == 0:
-        print "mkRamp(SIZE, DIRECTION, SLOPE, INTERCEPT, ORIGIN)"
-        print "first argument is required"
+        print("mkRamp(SIZE, DIRECTION, SLOPE, INTERCEPT, ORIGIN)")
+        print("first argument is required")
         exit(1)
     else:
         sz = args[0]
         if isinstance(sz, (int)):
             sz = (sz, sz)
         elif not isinstance(sz, (tuple)):
-            print "first argument must be a two element tuple or an integer"
+            print("first argument must be a two element tuple or an integer")
             exit(1)
 
     # OPTIONAL args:
@@ -697,8 +697,8 @@ def mkRamp(*args):
     xinc = slope * math.cos(direction)
     yinc = slope * math.sin(direction)
 
-    [xramp, yramp] = numpy.meshgrid( xinc * (numpy.array(range(sz[1]))-origin[1]),
-                                  yinc * (numpy.array(range(sz[0]))-origin[0]) )
+    [xramp, yramp] = numpy.meshgrid( xinc * (numpy.array(list(range(sz[1])))-origin[1]),
+                                  yinc * (numpy.array(list(range(sz[0])))-origin[0]) )
 
     res = intercept + xramp + yramp
 
@@ -1568,8 +1568,8 @@ def nextSz(size, sizeList):
 def mkImpulse(*args):
     # create an image that is all zeros except for an impulse
     if(len(args) == 0):
-        print "mkImpulse(size, origin, amplitude)"
-        print "first input parameter is required"
+        print("mkImpulse(size, origin, amplitude)")
+        print("first input parameter is required")
         return
     
     if(isinstance(args[0], int)):
@@ -1577,7 +1577,7 @@ def mkImpulse(*args):
     elif(isinstance(args[0], tuple)):
         sz = args[0]
     else:
-        print "size parameter must be either an integer or a tuple"
+        print("size parameter must be either an integer or a tuple")
         return
 
     if(len(args) > 1):
@@ -1606,7 +1606,7 @@ def mkImpulse(*args):
 def steer2HarmMtx(*args):
 
     if len(args) == 0:
-        print "Error: first parameter 'harmonics' is required."
+        print("Error: first parameter 'harmonics' is required.")
         return
     
     if len(args) > 0:
@@ -1617,19 +1617,19 @@ def steer2HarmMtx(*args):
     if len(args) > 1:
         angles = args[1]
     else:
-        angles = numpy.pi * numpy.array(range(numh)) / numh
+        angles = numpy.pi * numpy.array(list(range(numh))) / numh
         
     if len(args) > 2:
-        if isinstance(args[2], basestring):
+        if isinstance(args[2], str):
             if args[2] == 'even' or args[2] == 'EVEN':
                 evenorodd = 0
             elif args[2] == 'odd' or args[2] == 'ODD':
                 evenorodd = 1
             else:
-                print "Error: only 'even' and 'odd' are valid entries for the third input parameter."
+                print("Error: only 'even' and 'odd' are valid entries for the third input parameter.")
                 return
         else:
-            print "Error: third input parameter must be a string (even/odd)."
+            print("Error: third input parameter must be a string (even/odd).")
     else:
         evenorodd = 0
 
@@ -1653,7 +1653,7 @@ def steer2HarmMtx(*args):
 
     r = numpy.rank(imtx)
     if r != numh and r != angles.shape[0]:
-        print "Warning: matrix is not full rank"
+        print("Warning: matrix is not full rank")
 
     mtx = numpy.linalg.pinv(imtx)
     
@@ -1692,7 +1692,7 @@ def rcosFn(*args):
 
     sz = 256   # arbitrary!
 
-    X = numpy.pi * numpy.array(range(-sz-1,2)) / (2*sz)
+    X = numpy.pi * numpy.array(list(range(-sz-1,2))) / (2*sz)
 
     Y = values[0] + (values[1]-values[0]) * numpy.cos(X)**2;
 
@@ -1715,8 +1715,8 @@ def mkAngle(*args):
         if not isinstance(sz, tuple):
             sz = (sz, sz)
     else:
-        print "Error: first input parameter 'size' is required!"
-        print "makeAngle(size, phase, origin)"
+        print("Error: first input parameter 'size' is required!")
+        print("makeAngle(size, phase, origin)")
         return
 
     # ------------------------------------------------------------
@@ -1734,8 +1734,8 @@ def mkAngle(*args):
 
     #------------------------------------------------------------------
 
-    (xramp, yramp) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-origin[1], 
-                                 (numpy.array(range(1,sz[0]+1)))-origin[0])
+    (xramp, yramp) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-origin[1], 
+                                 (numpy.array(list(range(1,sz[0]+1))))-origin[0])
     xramp = numpy.array(xramp)
     yramp = numpy.array(yramp)
 
@@ -1754,7 +1754,7 @@ def mkAngle(*args):
 def modulateFlip(*args):
 
     if len(args) == 0:
-        print "Error: filter input parameter required."
+        print("Error: filter input parameter required.")
         return
 
     lfilt = args[0]
@@ -1764,13 +1764,13 @@ def modulateFlip(*args):
     elif lfilt.shape[0] == 1:
         lfilt = lfilt.reshape(lfilt.shape[1], 1)
     elif len(lfilt.shape) > 2 or lfilt.shape[1] != 1:
-        print 'Error: only 1D input supported.'
+        print('Error: only 1D input supported.')
         return
 
     sz = len(lfilt)
     sz2 = numpy.ceil(sz/2.0);
 
-    ind = numpy.array(range(sz-1,-1,-1))
+    ind = numpy.array(list(range(sz-1,-1,-1)))
 
     hfilt = lfilt[ind].T * (-1)**((ind+1)-sz2)
 
@@ -1794,7 +1794,7 @@ def modulateFlip(*args):
 # function res = blurDn(im, nlevs, filt)
 def blurDn(*args):
     if len(args) == 0:
-        print "Error: image input parameter required."
+        print("Error: image input parameter required.")
         return
 
     im = numpy.array(args[0])
@@ -1807,7 +1807,7 @@ def blurDn(*args):
 
     if len(args) > 2:
         filt = args[2]
-        if isinstance(filt, basestring):
+        if isinstance(filt, str):
             filt = namedFilter(filt)
     else:
         filt = namedFilter('binom5')
@@ -1827,7 +1827,7 @@ def blurDn(*args):
             # 1D image
             if len(filt.shape) > 1 and (filt.shape[1]!=1 and filt.shape[2]!=1):
                 # >1D filter
-                print 'Error: Cannot apply 2D filter to 1D signal'
+                print('Error: Cannot apply 2D filter to 1D signal')
                 return
             # orient filter and image correctly
             if im.shape[0] == 1:
@@ -1869,8 +1869,8 @@ def blur(*args):
 
     # REQUIRED ARG:
     if len(args) == 0:
-        print "blur(IM, LEVELS, FILT)"
-        print "first argument is required"
+        print("blur(IM, LEVELS, FILT)")
+        print("first argument is required")
         exit(1)
     else:
         im = numpy.array(args[0])
@@ -1882,7 +1882,7 @@ def blur(*args):
         nlevs = 1
 
     if len(args) > 2:
-        if isinstance(args[2], basestring):
+        if isinstance(args[2], str):
             filt = namedFilter(args[2])
         else:
             filt = numpy.array(args[2])
@@ -1900,7 +1900,7 @@ def blur(*args):
         if len(im.shape) == 1 or im.shape[0] == 1 or im.shape[1] == 1: 
             # 1D image
             if len(filt) == 2 and (numpy.asarray(filt.shape) != 1).any():
-                print 'Error: can not apply 2D filter to 1D signal'
+                print('Error: can not apply 2D filter to 1D signal')
                 return
             
             imIn = corrDn(im, filt, 'reflect1', len(im))
@@ -1941,7 +1941,7 @@ def rconv2(*args):
     #       1   (DIM/2)+1  
     
     if len(args) < 2:
-        print "Error: two matrices required as input parameters"
+        print("Error: two matrices required as input parameters")
         return
 
     if len(args) == 2:
@@ -1956,7 +1956,7 @@ def rconv2(*args):
         large = args[1]
         small = args[0]
     else:
-        print 'one arg must be larger than the other in both dimensions!'
+        print('one arg must be larger than the other in both dimensions!')
         return
 
     ly = large.shape[0]
@@ -1992,7 +1992,7 @@ def rconv2(*args):
 # compute minimum and maximum values of input matrix, returning them as tuple
 def range2(*args):
     if not numpy.isreal(args[0]).all():
-        print 'Error: matrix must be real-valued'
+        print('Error: matrix must be real-valued')
 
     return (args[0].min(), args[0].max())
 
@@ -2017,7 +2017,7 @@ def var2(*args):
 #  MEAN (optional) and VAR (optional) make the computation faster.
 def kurt2(*args):
     if len(args) == 0:
-        print 'Error: input matrix is required'
+        print('Error: input matrix is required')
 
     if len(args) < 2:
         mn = args[0].mean()
@@ -2045,16 +2045,16 @@ def kurt2(*args):
 def imStats(*args):
 
     if len(args) == 0:
-        print 'Error: at least one input image is required'
+        print('Error: at least one input image is required')
         return
     elif len(args) == 1 and not numpy.isreal(args[0]).all():
-        print 'Error: input images must be real-valued matrices'
+        print('Error: input images must be real-valued matrices')
         return
     elif len(args) == 2 and ( not numpy.isreal(args[0]).all() or not numpy.isreal(args[1]).all()):
-        print 'Error: input images must be real-valued matrices'
+        print('Error: input images must be real-valued matrices')
         return
     elif len(args) > 2:
-        print 'Error: maximum of two input images allowed'
+        print('Error: maximum of two input images allowed')
         return
 
     if len(args) == 2:
@@ -2066,18 +2066,18 @@ def imStats(*args):
             snr = numpy.inf
         else:
             snr = 10 * numpy.log10(var2(args[0])/v)
-        print 'Difference statistics:'
-        print '  Range: [%d, %d]' % (mn, mx)
-        print '  Mean: %f,  Stdev (rmse): %f,  SNR (dB): %f' % (mean, numpy.sqrt(v), snr)
+        print('Difference statistics:')
+        print('  Range: [%d, %d]' % (mn, mx))
+        print('  Mean: %f,  Stdev (rmse): %f,  SNR (dB): %f' % (mean, numpy.sqrt(v), snr))
     else:
         (mn, mx) = range2(args[0])
         mean = args[0].mean()
         var = var2(args[0])
         stdev = numpy.sqrt(var.real) + numpy.sqrt(var.imag)
         kurt = kurt2(args[0], mean, stdev**2)
-        print 'Image statistics:'
-        print '  Range: [%f, %f]' % (mn, mx)
-        print '  Mean: %f,  Stdev: %f,  Kurtosis: %f' % (mean, stdev, kurt)
+        print('Image statistics:')
+        print('  Range: [%f, %f]' % (mn, mx))
+        print('  Mean: %f,  Stdev: %f,  Kurtosis: %f' % (mean, stdev, kurt))
         
 # makes image the same as read in by matlab
 def correctImage(img):
@@ -2114,12 +2114,12 @@ def shift(mtx, offset):
 # Eero Simoncelli, 6/96.  Ported to Python by Rob Young, 5/14.
 def mkR(*args):
     if len(args) == 0:
-        print 'Error: first input parameter is required!'
+        print('Error: first input parameter is required!')
         return
     else:
         sz = args[0]
 
-    if isinstance(sz, (int, long)) or len(sz) == 1:
+    if isinstance(sz, int) or len(sz) == 1:
         sz = (sz, sz)
 
     # -----------------------------------------------------------------
@@ -2137,8 +2137,8 @@ def mkR(*args):
 
     # -----------------------------------------------------------------
 
-    (xramp2, yramp2) = numpy.meshgrid(numpy.array(range(1,sz[1]+1))-origin[1], 
-                                   numpy.array(range(1,sz[0]+1))-origin[0])
+    (xramp2, yramp2) = numpy.meshgrid(numpy.array(list(range(1,sz[1]+1)))-origin[1], 
+                                   numpy.array(list(range(1,sz[0]+1)))-origin[0])
 
     res = (xramp2**2 + yramp2**2)**(expt/2.0)
     
@@ -2155,9 +2155,9 @@ def mkR(*args):
 
 def mkFract(*args):
     if len(args) == 0:
-        print 'Error: input parameter dims required'
+        print('Error: input parameter dims required')
     else:
-        if isinstance(args[0], (int, long)) or len(args[0]) == 1:
+        if isinstance(args[0], int) or len(args[0]) == 1:
             dims = (args[0], args[0])
         elif args[0] == 1:
             dims = (args[1], args[1])
@@ -2185,7 +2185,7 @@ def mkFract(*args):
 
     #if any(max(max(abs(fres.imag))) > 1e-10):
     if abs(fres.imag).max() > 1e-10:
-        print 'Symmetry error in creating fractal'
+        print('Symmetry error in creating fractal')
     else:
         res = numpy.real(fres)
         res = res / numpy.sqrt(var2(res))
@@ -2213,7 +2213,7 @@ def mkFract(*args):
 def steer(*args):
     
     if len(args) < 2:
-        print 'Error: input parameters basis and angle are required!'
+        print('Error: input parameters basis and angle are required!')
         return
 
     basis = args[0]
@@ -2222,34 +2222,34 @@ def steer(*args):
 
     #if ( any(size(angle) ~= [size(basis,1) 1]) & any(size(angle) ~= [1 1]) )
     angle = args[1]
-    if isinstance(angle, (int, long, float)):
+    if isinstance(angle, (int, float)):
         angle = numpy.array([angle])
     else:
         if angle.shape[0] != basis.shape[0] or angle.shape[1] != 1:
-            print 'ANGLE must be a scalar, or a column vector the size of the basis elements'
+            print('ANGLE must be a scalar, or a column vector the size of the basis elements')
             return
 
     # If HARMONICS are not passed, assume derivatives.
     if len(args) < 3:
         if num%2 == 0:
-            harmonics = numpy.array(range(num/2))*2+1
+            harmonics = numpy.array(list(range(num/2)))*2+1
         else:
-            harmonics = numpy.array(range((15+1)/2))*2
+            harmonics = numpy.array(list(range((15+1)/2)))*2
     else:
         harmonics = args[2]
     if len(harmonics.shape) == 1 or harmonics.shape[0] == 1:
         harmonics = harmonics.reshape(harmonics.shape[0], 1)
     elif harmonics.shape[0] != 1 and harmonics.shape[1] != 1:
-        print 'Error: input parameter HARMONICS must be 1D!'
+        print('Error: input parameter HARMONICS must be 1D!')
         return
 
     if 2*harmonics.shape[0] - (harmonics == 0).sum() != num:
-        print 'harmonics list is incompatible with basis size!'
+        print('harmonics list is incompatible with basis size!')
         return
 
     # If STEERMTX not passed, assume evenly distributed cosine-phase filters:
     if len(args) < 4:
-        steermtx = steer2HarmMtx(harmonics, numpy.pi*numpy.array(range(num))/num, 
+        steermtx = steer2HarmMtx(harmonics, numpy.pi*numpy.array(list(range(num)))/num, 
                                  'even')
     else:
         steermtx = args[3]
@@ -2257,12 +2257,12 @@ def steer(*args):
     steervect = numpy.zeros((angle.shape[0], num))
     arg = angle * harmonics[numpy.nonzero(harmonics)[0]].T
     if all(harmonics):
-	steervect[:, range(0,num,2)] = numpy.cos(arg)
-	steervect[:, range(1,num,2)] = numpy.sin(arg)
+	steervect[:, list(range(0,num,2))] = numpy.cos(arg)
+	steervect[:, list(range(1,num,2))] = numpy.sin(arg)
     else:
 	steervect[:, 1] = numpy.ones((arg.shape[0],1))
-	steervect[:, range(0,num,2)] = numpy.cos(arg)
-	steervect[:, range(1,num,2)] = numpy.sin(arg)
+	steervect[:, list(range(0,num,2))] = numpy.cos(arg)
+	steervect[:, list(range(1,num,2))] = numpy.sin(arg)
 
     steervect = numpy.dot(steervect,steermtx)
 
@@ -2277,33 +2277,33 @@ def steer(*args):
 def showIm_old(*args):
     # check and set input parameters
     if len(args) == 0:
-        print "showIm( matrix, range, zoom, label, nshades )"
-        print "  matrix is string. It should be the name of a 2D array."
-        print "  range is a two element tuple.  It specifies the values that "
-        print "    map to the min and max colormap values.  Passing a value "
-        print "    of 'auto' (default) sets range=[min,max].  'auto2' sets "
-        print "    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets "
-        print "    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th "
-        print "    percientile value of the sorted matix samples, and p2 is "
-        print "    the 90th percentile value."
-        print "  zoom specifies the number of matrix samples per screen pixel."
-        print "    It will be rounded to an integer, or 1 divided by an "
-        print "    integer."
+        print("showIm( matrix, range, zoom, label, nshades )")
+        print("  matrix is string. It should be the name of a 2D array.")
+        print("  range is a two element tuple.  It specifies the values that ")
+        print("    map to the min and max colormap values.  Passing a value ")
+        print("    of 'auto' (default) sets range=[min,max].  'auto2' sets ")
+        print("    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets ")
+        print("    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th ")
+        print("    percientile value of the sorted matix samples, and p2 is ")
+        print("    the 90th percentile value.")
+        print("  zoom specifies the number of matrix samples per screen pixel.")
+        print("    It will be rounded to an integer, or 1 divided by an ")
+        print("    integer.")
         #print "    A value of 'same' or 'auto' (default) causes the "
         #print "    zoom value to be chosen automatically to fit the image into"
         #print "    the current axes."
         #print "    A value of 'full' fills the axis region "
         #print "    (leaving no room for labels)."
-        print "  label - A string that is used as a figure title."
-        print "  NSHADES (optional) specifies the number of gray shades, "
-        print "    and defaults to the size of the current colormap. "
+        print("  label - A string that is used as a figure title.")
+        print("  NSHADES (optional) specifies the number of gray shades, ")
+        print("    and defaults to the size of the current colormap. ")
 
     if len(args) > 0:   # matrix entered
         matrix = numpy.array(args[0])
     #print 'showIm range %f %f' % (matrix.min(), matrix.max())
 
     if len(args) > 1:   # range entered
-        if isinstance(args[1], basestring):
+        if isinstance(args[1], str):
             if args[1] is "auto":
                 imRange = ( numpy.amin(matrix), numpy.amax(matrix) )
             elif args[1] is "auto2":
@@ -2316,10 +2316,10 @@ def showIm_old(*args):
                 p2 = scipy.stats.scoreatpercentile(numpy.hstack(matrix), 90)
                 imRange = (p1-(p2-p1)/8.0, p2+(p2-p1)/8.0)
             else:
-                print "Error: range of %s is not recognized." % args[1]
-                print "       please use a two element tuple or "
-                print "       'auto', 'auto2' or 'auto3'"
-                print "       enter 'showIm' for more info about options"
+                print("Error: range of %s is not recognized." % args[1])
+                print("       please use a two element tuple or ")
+                print("       'auto', 'auto2' or 'auto3'")
+                print("       enter 'showIm' for more info about options")
                 return
         else:
             imRange = args[1][0], args[1][1]
@@ -2363,7 +2363,7 @@ def showIm_old(*args):
     
     # make colormap
     incr = (256/nshades)+1
-    colors = range(0,255,(256/nshades)+1)
+    colors = list(range(0,255,(256/nshades)+1))
     colors[-1] = 255
     colctr = -1
     for i in range(256):
@@ -2375,9 +2375,9 @@ def showIm_old(*args):
     # zoom
     #dims = (matrix.shape[0]*zoom, matrix.shape[1]*zoom)
     dims = (nRows*zoom, nCols*zoom)
-    print 'dims'
-    print dims
-    print 'nRows=%d nCols=%d' % (nRows, nCols)
+    print('dims')
+    print(dims)
+    print('nRows=%d nCols=%d' % (nRows, nCols))
     #qim = qim.scaled(dims[0], dims[1])
     qim = qim.scaled(nCols, nRows)
     #pixmap = QtGui.QPixmap()
@@ -2417,33 +2417,33 @@ def showIm_old(*args):
 def showIm(*args):
     # check and set input parameters
     if len(args) == 0:
-        print "showIm( matrix, range, zoom, label, nshades )"
-        print "  matrix is string. It should be the name of a 2D array."
-        print "  range is a two element tuple.  It specifies the values that "
-        print "    map to the min and max colormap values.  Passing a value "
-        print "    of 'auto' (default) sets range=[min,max].  'auto2' sets "
-        print "    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets "
-        print "    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th "
-        print "    percientile value of the sorted matix samples, and p2 is "
-        print "    the 90th percentile value."
-        print "  zoom specifies the number of matrix samples per screen pixel."
-        print "    It will be rounded to an integer, or 1 divided by an "
-        print "    integer."
+        print("showIm( matrix, range, zoom, label, nshades )")
+        print("  matrix is string. It should be the name of a 2D array.")
+        print("  range is a two element tuple.  It specifies the values that ")
+        print("    map to the min and max colormap values.  Passing a value ")
+        print("    of 'auto' (default) sets range=[min,max].  'auto2' sets ")
+        print("    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets ")
+        print("    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th ")
+        print("    percientile value of the sorted matix samples, and p2 is ")
+        print("    the 90th percentile value.")
+        print("  zoom specifies the number of matrix samples per screen pixel.")
+        print("    It will be rounded to an integer, or 1 divided by an ")
+        print("    integer.")
         #print "    A value of 'same' or 'auto' (default) causes the "
         #print "    zoom value to be chosen automatically to fit the image into"
         #print "    the current axes."
         #print "    A value of 'full' fills the axis region "
         #print "    (leaving no room for labels)."
-        print "  label - A string that is used as a figure title."
-        print "  NSHADES (optional) specifies the number of gray shades, "
-        print "    and defaults to the size of the current colormap. "
+        print("  label - A string that is used as a figure title.")
+        print("  NSHADES (optional) specifies the number of gray shades, ")
+        print("    and defaults to the size of the current colormap. ")
 
     if len(args) > 0:   # matrix entered
         matrix = numpy.array(args[0])
     #print 'showIm range %f %f' % (matrix.min(), matrix.max())
 
     if len(args) > 1:   # range entered
-        if isinstance(args[1], basestring):
+        if isinstance(args[1], str):
             if args[1] is "auto":
                 imRange = ( numpy.amin(matrix), numpy.amax(matrix) )
             elif args[1] is "auto2":
@@ -2456,10 +2456,10 @@ def showIm(*args):
                 p2 = scipy.stats.scoreatpercentile(numpy.hstack(matrix), 90)
                 imRange = (p1-(p2-p1)/8.0, p2+(p2-p1)/8.0)
             else:
-                print "Error: range of %s is not recognized." % args[1]
-                print "       please use a two element tuple or "
-                print "       'auto', 'auto2' or 'auto3'"
-                print "       enter 'showIm' for more info about options"
+                print("Error: range of %s is not recognized." % args[1])
+                print("       please use a two element tuple or ")
+                print("       'auto', 'auto2' or 'auto3'")
+                print("       enter 'showIm' for more info about options")
                 return
         else:
             imRange = args[1][0], args[1][1]
@@ -2483,24 +2483,24 @@ def showIm(*args):
 
     # create window
     #master = Tkinter.Tk()
-    master = Tkinter.Toplevel()
+    master = tkinter.Toplevel()
     master.title('showIm')
     canvas_width = matrix.shape[0] * zoom
     canvas_height = matrix.shape[1] * zoom
     master.geometry(str(canvas_width+20) + "x" + str(canvas_height+60) +
                     "+200+200")
     # put in top spacer
-    spacer = Tkinter.Label(master, text='').pack()
+    spacer = tkinter.Label(master, text='').pack()
     
     # create canvas
-    canvas = Tkinter.Canvas(master, width=canvas_width, height=canvas_height)
+    canvas = tkinter.Canvas(master, width=canvas_width, height=canvas_height)
     canvas.pack()
     #img = Image.fromarray(matrix)
     # FIX: shift matrix to 0.0-1.0 then to 0-255
     if (matrix < 0).any():
         matrix = matrix + math.fabs(matrix.min())
     matrix = (matrix / matrix.max()) * 255.0
-    print matrix.astype('uint8')[0,:]
+    print(matrix.astype('uint8')[0,:])
     img = PIL.Image.fromarray(matrix.astype('uint8'))
 
     # make colormap - works without range
@@ -2523,7 +2523,7 @@ def showIm(*args):
     #incr = int(numpy.ceil(float(imRange[1]-imRange[0]+1) / float(nshades)))
     incr = int(numpy.ceil(float(matrix.max()-matrix.min()+1) / float(nshades)))
     #colors = range(int(imRange[0]), int(imRange[1])+1, incr)
-    colors = range(int(matrix.min()), int(matrix.max())+1, incr)
+    colors = list(range(int(matrix.min()), int(matrix.max())+1, incr))
     colors[0] = 0
     colors[-1] = 255
     colctr = -1
@@ -2543,21 +2543,21 @@ def showIm(*args):
 
     # apply image to canvas
     imgPI = ImageTk.PhotoImage(img)    
-    canvas.create_image(0,0, anchor=Tkinter.NW, image=imgPI)
+    canvas.create_image(0,0, anchor=tkinter.NW, image=imgPI)
 
     # add labels
     rangeStr = 'Range: [%.1f, %.1f]' % (imRange[0], imRange[1])
-    rangeLabel = Tkinter.Label(master, text=rangeStr).pack()
+    rangeLabel = tkinter.Label(master, text=rangeStr).pack()
     dimsStr = 'Dims: [%d, %d] / %d' % (matrix.shape[0], matrix.shape[1], zoom)
-    dimsLabel = Tkinter.Label(master, text=dimsStr).pack()
+    dimsLabel = tkinter.Label(master, text=dimsStr).pack()
     
-    Tkinter.mainloop()
+    tkinter.mainloop()
 
 def corrDn(image = None, filt = None, edges = 'reflect1', step = (1,1), 
            start = (0,0), stop = None, result = None):
 
     if image == None or filt == None:
-        print 'Error: image and filter are required input parameters!'
+        print('Error: image and filter are required input parameters!')
         return
     else:
         image = image.copy()
@@ -2570,8 +2570,8 @@ def corrDn(image = None, filt = None, edges = 'reflect1', step = (1,1),
         stop = (image.shape[0], image.shape[1])
 
     if result == None:
-        rxsz = len(range(start[0], stop[0], step[0]))
-        rysz = len(range(start[1], stop[1], step[1]))
+        rxsz = len(list(range(start[0], stop[0], step[0])))
+        rysz = len(list(range(start[1], stop[1], step[1])))
         result = numpy.zeros((rxsz, rysz))
         
     if edges == 'circular':
@@ -2600,7 +2600,7 @@ def upConv(image = None, filt = None, edges = 'reflect1', step = (1,1),
            start = (0,0), stop = None, result = None):
 
     if image == None or filt == None:
-        print 'Error: image and filter are required input parameters!'
+        print('Error: image and filter are required input parameters!')
         return
     else:
         image = image.copy()
@@ -2619,7 +2619,7 @@ def upConv(image = None, filt = None, edges = 'reflect1', step = (1,1),
             filt = numpy.append(filt,0.0);
             filt = numpy.reshape(filt, (1, len(filt)))
         else:
-            print 'Even sized 2D filters not yet supported by upConv.'
+            print('Even sized 2D filters not yet supported by upConv.')
             return
 
     if stop == None and result == None:
@@ -2682,9 +2682,9 @@ def cconv2(*args):
     # Eero Simoncelli, 6/96.  Modified 2/97.  Python port by Rob Young, 8/15
     
     if len(args) < 2:
-        print 'Error: cconv2 requires two input matrices!'
-        print 'Usage: cconv2(matrix1, matrix2, center)'
-        print 'where center parameter is optional'
+        print('Error: cconv2 requires two input matrices!')
+        print('Usage: cconv2(matrix1, matrix2, center)')
+        print('where center parameter is optional')
         return
     else:
         a = numpy.array(args[0])
@@ -2702,7 +2702,7 @@ def cconv2(*args):
         large = b
         small = a
     else:
-        print 'Error: one matrix must be larger than the other in both dimensions!'
+        print('Error: one matrix must be larger than the other in both dimensions!')
         return
     
     ly = large.shape[0]
@@ -2748,8 +2748,8 @@ def clip(*args):
     # ported to Python by Rob Young, 8/15
     
     if len(args) == 0 or len(args) > 3:
-        print 'Usage: clip(im, minVal or Range, maxVal)'
-        print 'first input parameter is required'
+        print('Usage: clip(im, minVal or Range, maxVal)')
+        print('first input parameter is required')
         return
         
     im = numpy.array(args[0])
@@ -2758,7 +2758,7 @@ def clip(*args):
         minVal = 0;
         maxVal = 1;
     elif len(args) == 2:
-        if isinstance(args[1], (int, long, float)):
+        if isinstance(args[1], (int, float)):
             minVal = args[1]
             maxVal = args[1]+1
         else:
@@ -2769,7 +2769,7 @@ def clip(*args):
         maxVal = args[2]
         
     if maxVal < minVal:
-        print 'Error: maxVal cannot be less than minVal!'
+        print('Error: maxVal cannot be less than minVal!')
         return
 
     im[numpy.where(im < minVal)] = minVal
@@ -2781,7 +2781,7 @@ def clip(*args):
 # used in histo so we can unit test against matlab code
 # numpy version rounds to closest even number to remove bias
 def round(arr):
-    if isinstance(arr, (int, float, long)):
+    if isinstance(arr, (int, float)):
         arr = roundVal(arr)
     else:
         for i in range(len(arr)):
@@ -2822,8 +2822,8 @@ def histo(*args):
     #       called histo.c
 
     if len(args) == 0 or len(args) > 3:
-        print 'Usage: histo(mtx, nbins, binCtr)'
-        print 'first argument is required'
+        print('Usage: histo(mtx, nbins, binCtr)')
+        print('first argument is required')
         return
     else:
         mtx = args[0]
@@ -2844,7 +2844,7 @@ def histo(*args):
             tmpNbins = ( round(float(mx-binCtr) / float(binSize)) - 
                          round(float(mn-binCtr) / float(binSize)) )
             if tmpNbins != args[1]:
-                print 'Warning: Using %d bins instead of requested number (%d)' % (tmpNbins, args[1])
+                print('Warning: Using %d bins instead of requested number (%d)' % (tmpNbins, args[1]))
     else:
         binSize = float(mx-mn) / 101.0
 
@@ -2857,7 +2857,7 @@ def histo(*args):
     # numpy.histogram uses bin edges, not centers like Matlab's hist
     #bins = firstBin + binSize * numpy.array(range(tmpNbins+1))
     # compute bin edges
-    binsE = firstEdge + binSize * numpy.array(range(tmpNbins+1))
+    binsE = firstEdge + binSize * numpy.array(list(range(tmpNbins+1)))
     
     [N, X] = numpy.histogram(mtx, binsE)
 
@@ -2911,7 +2911,7 @@ def factorial(*args):
     # EPS, 11/02, Python port by Rob Young, 10/15
 
     # if scalar input make it a single element array
-    if isinstance(args[0], (int, long, float)):
+    if isinstance(args[0], (int, float)):
         num = numpy.array([args[0]])
     else:
         num = numpy.array(args[0])
@@ -2995,8 +2995,8 @@ def imGradient(*args):
     # Python port by Rob Young, 10/15
     
     if len(args) == 0 or len(args) > 2:
-        print 'Usage: imGradient(image, edges)'
-        print "'edges' argument is optional"
+        print('Usage: imGradient(image, edges)')
+        print("'edges' argument is optional")
     elif len(args) == 1:
         edges = "dont-compute"
     elif len(args) == 2:
@@ -3019,8 +3019,8 @@ def skew2(*args):
     #  MEAN (optional) and VAR (optional) make the computation faster.
 
     if len(args) == 0:
-        print 'Usage: skew2(matrix, mean, variance)'
-        print 'mean and variance arguments are optional'
+        print('Usage: skew2(matrix, mean, variance)')
+        print('mean and variance arguments are optional')
     else:
         mtx = numpy.array(args[0])
 
@@ -3060,8 +3060,8 @@ def upBlur(*args):
     # REQUIRED ARGS
     
     if len(args) == 0:
-        print 'Usage: upBlur(image, levels, filter)'
-        print 'first argument is required'
+        print('Usage: upBlur(image, levels, filter)')
+        print('first argument is required')
     else:
         im = numpy.array(args[0])
 
@@ -3080,7 +3080,7 @@ def upBlur(*args):
 
     #------------------------------------------------------------------
 
-    if isinstance(filt, basestring):
+    if isinstance(filt, str):
         filt = namedFilter(filt)
 
     if nlevs > 1:
@@ -3125,8 +3125,8 @@ def zconv2(*args):
     #----------------------------------------------------------------
     
     if len(args) < 2 or len(args) > 3:
-        print 'Usage: zconv2(matrix1, matrix2, center)'
-        print 'first two input parameters are required'
+        print('Usage: zconv2(matrix1, matrix2, center)')
+        print('first two input parameters are required')
     else:
         a = numpy.array(args[0])
         b = numpy.array(args[1])
@@ -3148,7 +3148,7 @@ def zconv2(*args):
         large = b
         small = a
     else:
-        print 'Error: one arg must be larger than the other in both dimensions!'
+        print('Error: one arg must be larger than the other in both dimensions!')
         return
         
     ly = large.shape[0]
diff --git a/pyramid.py b/pyramid.py
index 6fcdd3a..6f5c793 100644
--- a/pyramid.py
+++ b/pyramid.py
@@ -9,7 +9,7 @@ class pyramid:  # pyramid
 
     # constructor
     def __init__(self):
-        print "please specify type of pyramid to create (Gpry, Lpyr, etc.)"
+        print("please specify type of pyramid to create (Gpry, Lpyr, etc.)")
         return
 
     # methods
diff --git a/range2.py b/range2.py
index 078ae8e..ef86f73 100644
--- a/range2.py
+++ b/range2.py
@@ -5,6 +5,6 @@ def range2(*args):
         as tuple  '''
 
     if not numpy.isreal(args[0]).all():
-        print 'Error: matrix must be real-valued'
+        print('Error: matrix must be real-valued')
 
     return (args[0].min(), args[0].max())
diff --git a/rconv2.py b/rconv2.py
index e73902e..7a1bf99 100644
--- a/rconv2.py
+++ b/rconv2.py
@@ -13,7 +13,7 @@ def rconv2(*args):
               1   (DIM/2)+1  '''
     
     if len(args) < 2:
-        print "Error: two matrices required as input parameters"
+        print("Error: two matrices required as input parameters")
         return
 
     if len(args) == 2:
@@ -28,7 +28,7 @@ def rconv2(*args):
         large = args[1]
         small = args[0]
     else:
-        print 'one arg must be larger than the other in both dimensions!'
+        print('one arg must be larger than the other in both dimensions!')
         return
 
     ly = large.shape[0]
diff --git a/rcosFn.py b/rcosFn.py
index e40b96b..1b5355b 100644
--- a/rcosFn.py
+++ b/rcosFn.py
@@ -33,7 +33,7 @@ def rcosFn(*args):
 
     sz = 256   # arbitrary!
 
-    X = numpy.pi * numpy.array(range(-sz-1,2)) / (2*sz)
+    X = numpy.pi * numpy.array(list(range(-sz-1,2))) / (2*sz)
 
     Y = values[0] + (values[1]-values[0]) * numpy.cos(X)**2;
 
diff --git a/round.py b/round.py
index a672b91..393f355 100644
--- a/round.py
+++ b/round.py
@@ -1,12 +1,12 @@
 import numpy
-from roundVal import roundVal
+from .roundVal import roundVal
 
 def round(arr):
     ''' round equivalent to matlab function
         used in histo so we can unit test against matlab code
         numpy version rounds to closest even number to remove bias  '''
 
-    if isinstance(arr, (int, float, long)):
+    if isinstance(arr, (int, float)):
         arr = roundVal(arr)
     else:
         for i in range(len(arr)):
diff --git a/showIm.py b/showIm.py
index f178cfd..97bcf9a 100644
--- a/showIm.py
+++ b/showIm.py
@@ -1,40 +1,40 @@
 import numpy
-import ImageTk
+from PIL import ImageTk
 import PIL
 import scipy.stats
-import Tkinter
+import tkinter
 import math
-from round import round
+from .round import round
 
 def showIm(*args):
     # check and set input parameters
     if len(args) == 0:
-        print "showIm( matrix, range, zoom, label, nshades )"
-        print "  matrix is string. It should be the name of a 2D array."
-        print "  range is a two element tuple.  It specifies the values that "
-        print "    map to the min and max colormap values.  Passing a value "
-        print "    of 'auto' (default) sets range=[min,max].  'auto2' sets "
-        print "    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets "
-        print "    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th "
-        print "    percientile value of the sorted matix samples, and p2 is "
-        print "    the 90th percentile value."
-        print "  zoom specifies the number of matrix samples per screen pixel."
-        print "    It will be rounded to an integer, or 1 divided by an "
-        print "    integer."
+        print("showIm( matrix, range, zoom, label, nshades )")
+        print("  matrix is string. It should be the name of a 2D array.")
+        print("  range is a two element tuple.  It specifies the values that ")
+        print("    map to the min and max colormap values.  Passing a value ")
+        print("    of 'auto' (default) sets range=[min,max].  'auto2' sets ")
+        print("    range=[mean-2*stdev, mean+2*stdev].  'auto3' sets ")
+        print("    range=[p1-(p2-p1)/8, p2+(p2-p1)/8], where p1 is the 10th ")
+        print("    percientile value of the sorted matix samples, and p2 is ")
+        print("    the 90th percentile value.")
+        print("  zoom specifies the number of matrix samples per screen pixel.")
+        print("    It will be rounded to an integer, or 1 divided by an ")
+        print("    integer.")
         #print "    A value of 'same' or 'auto' (default) causes the "
         #print "    zoom value to be chosen automatically to fit the image into"
         #print "    the current axes."
         #print "    A value of 'full' fills the axis region "
         #print "    (leaving no room for labels)."
-        print "  label - A string that is used as a figure title."
-        print "  NSHADES (optional) specifies the number of gray shades, "
-        print "    and defaults to the size of the current colormap. "
+        print("  label - A string that is used as a figure title.")
+        print("  NSHADES (optional) specifies the number of gray shades, ")
+        print("    and defaults to the size of the current colormap. ")
 
     if len(args) > 0:   # matrix entered
         matrix = numpy.array(args[0])
 
     if len(args) > 1:   # range entered
-        if isinstance(args[1], basestring):
+        if isinstance(args[1], str):
             if args[1] is "auto":
                 imRange = ( numpy.amin(matrix), numpy.amax(matrix) )
             elif args[1] is "auto2":
@@ -47,10 +47,10 @@ def showIm(*args):
                 p2 = scipy.stats.scoreatpercentile(numpy.hstack(matrix), 90)
                 imRange = (p1-(p2-p1)/8.0, p2+(p2-p1)/8.0)
             else:
-                print "Error: range of %s is not recognized." % args[1]
-                print "       please use a two element tuple or "
-                print "       'auto', 'auto2' or 'auto3'"
-                print "       enter 'showIm' for more info about options"
+                print("Error: range of %s is not recognized." % args[1])
+                print("       please use a two element tuple or ")
+                print("       'auto', 'auto2' or 'auto3'")
+                print("       enter 'showIm' for more info about options")
                 return
         else:
             imRange = args[1][0], args[1][1]
@@ -74,29 +74,29 @@ def showIm(*args):
 
     # create window
     #master = Tkinter.Tk()
-    master = Tkinter.Toplevel()
+    master = tkinter.Toplevel()
     master.title('showIm')
     canvas_width = matrix.shape[0] * zoom
     canvas_height = matrix.shape[1] * zoom
     master.geometry(str(canvas_width+20) + "x" + str(canvas_height+60) +
                     "+200+200")
     # put in top spacer
-    spacer = Tkinter.Label(master, text='').pack()
+    spacer = tkinter.Label(master, text='').pack()
     
     # create canvas
-    canvas = Tkinter.Canvas(master, width=canvas_width, height=canvas_height)
+    canvas = tkinter.Canvas(master, width=canvas_width, height=canvas_height)
     canvas.pack()
     # shift matrix to 0.0-1.0 then to 0-255
     if (matrix < 0).any():
         matrix = matrix + math.fabs(matrix.min())
     matrix = (matrix / matrix.max()) * 255.0
-    print matrix.astype('uint8')[0,:]
+    print(matrix.astype('uint8')[0,:])
     img = PIL.Image.fromarray(matrix.astype('uint8'))
 
     # make colormap
     colorTable = [0] * 256
     incr = int(numpy.ceil(float(matrix.max()-matrix.min()+1) / float(nshades)))
-    colors = range(int(matrix.min()), int(matrix.max())+1, incr)
+    colors = list(range(int(matrix.min()), int(matrix.max())+1, incr))
     colors[0] = 0
     colors[-1] = 255
     colctr = -1
@@ -115,12 +115,12 @@ def showIm(*args):
 
     # apply image to canvas
     imgPI = ImageTk.PhotoImage(img)    
-    canvas.create_image(0,0, anchor=Tkinter.NW, image=imgPI)
+    canvas.create_image(0,0, anchor=tkinter.NW, image=imgPI)
 
     # add labels
     rangeStr = 'Range: [%.1f, %.1f]' % (imRange[0], imRange[1])
-    rangeLabel = Tkinter.Label(master, text=rangeStr).pack()
+    rangeLabel = tkinter.Label(master, text=rangeStr).pack()
     dimsStr = 'Dims: [%d, %d] / %d' % (matrix.shape[0], matrix.shape[1], zoom)
-    dimsLabel = Tkinter.Label(master, text=dimsStr).pack()
+    dimsLabel = tkinter.Label(master, text=dimsStr).pack()
     
-    Tkinter.mainloop()
+    tkinter.mainloop()
diff --git a/skew2.py b/skew2.py
index 98d1309..6b8050b 100644
--- a/skew2.py
+++ b/skew2.py
@@ -1,13 +1,13 @@
 import numpy
-from var2 import var2
+from .var2 import var2
 
 def skew2(*args):
     ''' Sample skew (third moment divided by variance^3/2) of a matrix.
         MEAN (optional) and VAR (optional) make the computation faster.  '''
 
     if len(args) == 0:
-        print 'Usage: skew2(matrix, mean, variance)'
-        print 'mean and variance arguments are optional'
+        print('Usage: skew2(matrix, mean, variance)')
+        print('mean and variance arguments are optional')
     else:
         mtx = numpy.array(args[0])
 
diff --git a/steer.py b/steer.py
index f50adc7..ef84483 100644
--- a/steer.py
+++ b/steer.py
@@ -1,5 +1,5 @@
 import numpy
-from steer2HarmMtx import steer2HarmMtx
+from .steer2HarmMtx import steer2HarmMtx
 
 def steer(*args):
     ''' Steer BASIS to the specfied ANGLE.  
@@ -22,7 +22,7 @@ def steer(*args):
         Eero Simoncelli, 7/96. Ported to Python by Rob Young, 5/14.  '''
     
     if len(args) < 2:
-        print 'Error: input parameters basis and angle are required!'
+        print('Error: input parameters basis and angle are required!')
         return
 
     basis = args[0]
@@ -30,19 +30,19 @@ def steer(*args):
     num = basis.shape[1]
 
     angle = args[1]
-    if isinstance(angle, (int, long, float)):
+    if isinstance(angle, (int, float)):
         angle = numpy.array([angle])
     else:
         if angle.shape[0] != basis.shape[0] or angle.shape[1] != 1:
-            print 'ANGLE must be a scalar, or a column vector the size of the basis elements'
+            print('ANGLE must be a scalar, or a column vector the size of the basis elements')
             return
 
     # If HARMONICS are not passed, assume derivatives.
     if len(args) < 3:
         if num%2 == 0:
-            harmonics = numpy.array(range(num/2))*2+1
+            harmonics = numpy.array(list(range(num/2)))*2+1
         else:
-            harmonics = numpy.array(range((15+1)/2))*2
+            harmonics = numpy.array(list(range((15+1)/2)))*2
     else:
         harmonics = args[2]
 
@@ -50,17 +50,17 @@ def steer(*args):
         # reshape to column matrix
         harmonics = harmonics.reshape(harmonics.shape[0], 1)
     elif harmonics.shape[0] != 1 and harmonics.shape[1] != 1:
-        print 'Error: input parameter HARMONICS must be 1D!'
+        print('Error: input parameter HARMONICS must be 1D!')
         return
 
     if 2*harmonics.shape[0] - (harmonics == 0).sum() != num:
-        print 'harmonics list is incompatible with basis size!'
+        print('harmonics list is incompatible with basis size!')
         return
 
     # If STEERMTX not passed, assume evenly distributed cosine-phase filters:
     if len(args) < 4:
         steermtx = steer2HarmMtx(harmonics,
-                                 numpy.pi*numpy.array(range(num))/num, 
+                                 numpy.pi*numpy.array(list(range(num)))/num, 
                                  'even')
     else:
         steermtx = args[3]
@@ -68,19 +68,19 @@ def steer(*args):
     steervect = numpy.zeros((angle.shape[0], num))
     arg = angle * harmonics[numpy.nonzero(harmonics)[0]].T
     if all(harmonics):
-	steervect[:, range(0,num,2)] = numpy.cos(arg)
-	steervect[:, range(1,num,2)] = numpy.sin(arg)
+        steervect[:, list(range(0,num,2))] = numpy.cos(arg)
+        steervect[:, list(range(1,num,2))] = numpy.sin(arg)
     else:
-	steervect[:, 1] = numpy.ones((arg.shape[0],1))
-	steervect[:, range(0,num,2)] = numpy.cos(arg)
-	steervect[:, range(1,num,2)] = numpy.sin(arg)
+        steervect[:, 1] = numpy.ones((arg.shape[0],1))
+        steervect[:, list(range(0,num,2))] = numpy.cos(arg)
+        steervect[:, list(range(1,num,2))] = numpy.sin(arg)
 
     steervect = numpy.dot(steervect,steermtx)
 
     if steervect.shape[0] > 1:
-	tmp = numpy.dot(basis, steervect)
-	res = sum(tmp).T
+        tmp = numpy.dot(basis, steervect)
+        res = sum(tmp).T
     else:
-	res = numpy.dot(basis, steervect.T)
+        res = numpy.dot(basis, steervect.T)
 
     return res
diff --git a/steer2HarmMtx.py b/steer2HarmMtx.py
index 48c731b..ba59116 100644
--- a/steer2HarmMtx.py
+++ b/steer2HarmMtx.py
@@ -11,7 +11,7 @@ def steer2HarmMtx(*args):
         mtx = steer2HarmMtx(harmonics, angles, evenorodd)  '''
 
     if len(args) == 0:
-        print "Error: first parameter 'harmonics' is required."
+        print("Error: first parameter 'harmonics' is required.")
         return
     
     if len(args) > 0:
@@ -22,19 +22,19 @@ def steer2HarmMtx(*args):
     if len(args) > 1:
         angles = args[1]
     else:
-        angles = numpy.pi * numpy.array(range(numh)) / numh
+        angles = numpy.pi * numpy.array(list(range(numh))) / numh
         
     if len(args) > 2:
-        if isinstance(args[2], basestring):
+        if isinstance(args[2], str):
             if args[2] == 'even' or args[2] == 'EVEN':
                 evenorodd = 0
             elif args[2] == 'odd' or args[2] == 'ODD':
                 evenorodd = 1
             else:
-                print "Error: only 'even' and 'odd' are valid entries for the third input parameter."
+                print("Error: only 'even' and 'odd' are valid entries for the third input parameter.")
                 return
         else:
-            print "Error: third input parameter must be a string (even/odd)."
+            print("Error: third input parameter must be a string (even/odd).")
     else:
         evenorodd = 0
 
@@ -58,7 +58,7 @@ def steer2HarmMtx(*args):
 
     r = numpy.rank(imtx)
     if r != numh and r != angles.shape[0]:
-        print "Warning: matrix is not full rank"
+        print("Warning: matrix is not full rank")
 
     mtx = numpy.linalg.pinv(imtx)
     
diff --git a/upBlur.py b/upBlur.py
index 7de84ff..f8bd39a 100644
--- a/upBlur.py
+++ b/upBlur.py
@@ -1,6 +1,6 @@
 import numpy
-from namedFilter import namedFilter
-from upConv import upConv
+from .namedFilter import namedFilter
+from .upConv import upConv
 
 def upBlur(*args):
     ''' RES = upBlur(IM, LEVELS, FILT)
@@ -20,8 +20,8 @@ def upBlur(*args):
     # REQUIRED ARGS
     
     if len(args) == 0:
-        print 'Usage: upBlur(image, levels, filter)'
-        print 'first argument is required'
+        print('Usage: upBlur(image, levels, filter)')
+        print('first argument is required')
     else:
         im = numpy.array(args[0])
 
@@ -40,7 +40,7 @@ def upBlur(*args):
 
     #------------------------------------------------------------------
 
-    if isinstance(filt, basestring):
+    if isinstance(filt, str):
         filt = namedFilter(filt)
 
     if nlevs > 1:
diff --git a/upConv.py b/upConv.py
index 52b36d9..16fc74b 100644
--- a/upConv.py
+++ b/upConv.py
@@ -10,7 +10,7 @@ def upConv(image = None, filt = None, edges = 'reflect1', step = (1,1),
            start = (0,0), stop = None, result = None):
 
     if image is None or filt is None:
-        print 'Error: image and filter are required input parameters!'
+        print('Error: image and filter are required input parameters!')
         return
     else:
         image = image.copy()
@@ -29,7 +29,7 @@ def upConv(image = None, filt = None, edges = 'reflect1', step = (1,1),
             filt = numpy.append(filt,0.0);
             filt = numpy.reshape(filt, (1, len(filt)))
         else:
-            print 'Even sized 2D filters not yet supported by upConv.'
+            print('Even sized 2D filters not yet supported by upConv.')
             return
 
     if stop is None and result is None:
diff --git a/wrap.c b/wrap.c
index 3651b6b..de23af6 100755
--- a/wrap.c
+++ b/wrap.c
@@ -15,7 +15,8 @@
 */
 
 #include <stdlib.h>
-
+#include <stdio.h>
+#include <math.h>
 #include "convolve.h"
 
 /*
@@ -274,7 +275,55 @@ int internal_wrap_expand(image, filt, x_fdim, y_fdim,
   return(0);
   } /* end of internal_wrap_expand */
 
+void internal_pointop (im, res, size, lut, lutsize, origin, increment, warnings)
+  register double *im, *res, *lut;
+  register double origin, increment;
+  register int size, lutsize, warnings;
+  {
+  register int i, index;
+  register double pos;
+  register int l_unwarned = warnings;
+  register int r_unwarned = warnings;
+
+  lutsize = lutsize - 2;	/* Maximum index value */
 
+  /* printf("size=%d origin=%f lutsize=%d increment=%f\n",size, origin, lutsize,
+     increment); */
+
+  if (increment > 0)
+    for (i=0; i<size; i++)
+      {
+	pos = (im[i] - origin) / increment;
+	index = (int) pos;   /* Floor */
+	if (index < 0)
+	  {
+	    index = 0;
+	    if (l_unwarned)
+	      {
+		printf("Warning: Extrapolating to left of lookup table...\n");
+		l_unwarned = 0;
+	      }
+	  }
+	else if (index > lutsize)
+	  {
+	    index = lutsize;
+	    if (r_unwarned)
+	      {
+		printf("Warning: Extrapolating to right of lookup table...\n");
+		r_unwarned = 0;
+	      }
+	  }
+	res[i] = lut[index] + (lut[index+1] - lut[index]) * (pos - index);
+	if(isnan(res[i]))
+	  printf("**NAN: lut[%d]=%f lut[%d]=%f pos=%f index=%d\n", index,
+		 lut[index], index+1, lut[index+1], pos, index);
+      }
+  else
+    for (i=0; i<size; i++){
+      res[i] = *lut;
+      /*printf("res[%d]=%f\n", i, res[i]);*/
+    }
+  }
 
 /* Local Variables: */
 /* buffer-read-only: t */
diff --git a/wrapConv.so b/wrapConv.so
new file mode 100644
index 0000000..f37ae76
Binary files /dev/null and b/wrapConv.so differ
diff --git a/zconv2.py b/zconv2.py
index f64ef09..6cf3c46 100644
--- a/zconv2.py
+++ b/zconv2.py
@@ -20,8 +20,8 @@ def zconv2(*args):
     #----------------------------------------------------------------
     
     if len(args) < 2 or len(args) > 3:
-        print 'Usage: zconv2(matrix1, matrix2, center)'
-        print 'first two input parameters are required'
+        print('Usage: zconv2(matrix1, matrix2, center)')
+        print('first two input parameters are required')
         return
     else:
         a = numpy.array(args[0])
@@ -44,7 +44,7 @@ def zconv2(*args):
         large = b
         small = a
     else:
-        print 'Error: one arg must be larger than the other in both dimensions!'
+        print('Error: one arg must be larger than the other in both dimensions!')
         return
         
     ly = large.shape[0]