new

2 files changed, 1154 insertions, 400 deletions

diff --git a/bin/3Drotate.py b/bin/3Drotate.py
index 7669ffa..a94a7a9 100755
--- a/bin/3Drotate.py
+++ b/bin/3Drotate.py
@@ -1,11 +1,12 @@
 #!/usr/bin/python2.7
-import getopt
+import re
 import numpy
 import sys
+import os
 from subprocess import Popen, PIPE
 
 #{{{ usage
-USAGE="""
+USAGE = """
 USAGE: 3Drotate option=value infile outfile
 USAGE: 3Drotate [-h or -help]
 
@@ -46,8 +47,9 @@ vp       value     virtual-pixel method; any valid IM virtual-pixel method;
 
 NAME: 3DROTATE
 
-PURPOSE: To apply a perspective distortion to an image by providing rotation angles,
-zoom, offsets, background color, perspective exaggeration and auto zoom/centering.
+PURPOSE: To apply a perspective distortion to an image by
+providing rotation angles, zoom, offsets, background color,
+perspective exaggeration and auto zoom/centering.
 
 DESCRIPTION: 3DROTATE applies a perspective distortion to an image
 by providing any combination of three optional rotation angle:
@@ -186,23 +188,27 @@ skycolor = "black"
 vp = "background"
 
 # set directory for temporary files
-dir = "."    # suggestions are dir="." or dir="/tmp"
+directory = "."    # suggestions are dir="." or dir="/tmp"
 
 pi = numpy.math.pi
 #}}}
 
+
 # function to do dot product of 2 three element vectors
 def dot(vec1, vec2):
     return numpy.dot(vec1, vec2)
 
-# function to do 3x3 matrix multiply M x N where input are rows of each matrix; M1 M2 M3 N1 N2 N3
+
+# function to do 3x3 matrix multiply M x N where input
+# are rows of each matrix; M1 M2 M3 N1 N2 N3
 def matmul3(m0, m1, m2, n0, n1, n2):
     return numpy.matmul(
       [m0, m1, m2], [n0, n1, n2]
     )
 
+
 # function to project points from input to output domain
-def forwardProjrect(mat, ii, jj):
+def forwardProject(mat, ii, jj):
     #takes in a 3x3 matrix and scalars ii and jj
     #returns two scalars, uu and vv
     numu = mat[0][0] * ii + mat[0][1] * jj + mat[0][2]
@@ -212,10 +218,17 @@ def forwardProjrect(mat, ii, jj):
     vv = numv / den
     return [uu, vv]
 
+
 def errMsg(s):
     sys.stderr.write("%s\n")
     sys.exit(1)
 
+
+def usage():
+    sys.stderr.write(USAGE)
+    sys.exit(1)
+
+
 def inverseProject(mat, uu, vv):
     #note this function is more exact in python version
     numi = (mat[0][0] * uu) + (mat[0][1] * vv) + mat[0][2]
@@ -226,453 +239,478 @@ def inverseProject(mat, uu, vv):
     #FIXME chop off decimal for it to be like bash version
     return [ii, jj]
 
+
 # function to invert a 3 x 3 matrix using method of adjoint
 # inverse is the transpose of the matrix of cofactors divided by the determinant
 def M3inverse(mat):
     return numpy.linalg.inv(mat)
 
+
 def call_cmd(s):
     cmd = s.split(" ")
-    ps = Popen(cmd, stdout=PIPE, stderr=PIPE)
+    process = Popen(cmd, stdout=PIPE, stderr=PIPE)
     out, err = process.communicate()
     errcode = process.returncode
     if errcode > 0:
         errMsg(err)
     return out
 
-#
+
 # get input image size
 def imagesize(filepath):
     width = int(call_cmd("identify -format %w %s" % (filepath)))
     height = int(call_cmd("identify -format %h %s" % (filepath)))
     return width, height
 
+
 # test for correct number of arguments and get values
 len_args = len(sys.argv)
+
 if len_args == 1:
-   usage2()
-   sys.exit(1)
+    usage()
 elif (len_args > 15):
     errMsg("--- TOO MANY ARGUMENTS WERE PROVIDED ---")
-else:
-    for arg in sys.argv:
-        if arg == "-h":
-            usage2()
+
+for arg in sys.argv:
+    if arg == "-h":
+        usage()
+    if arg == "-":
+        break  # FIXME
+    test = re.match(r'pan=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        pan = float(test.group(1))
+        if (pan > 180) or (pan < -180):
+            errMsg(
+                "PAN=%s MUST BE GREATER THAN -180 AND LESS THAN +180"
+            ) % pan
+            sys.exit(1)
+        panang = numpy.multiply(numpy.pi, numpy.divide(pan, 180))
+        sinpan = numpy.sin(panang)
+        cospan = numpy.cos(panang)
+        Rp0 = [cospan, 0, sinpan]
+        Rp1 = [0, 1, 0]
+        Rp2 = [-sinpan, 0, cospan]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rp0, Rp1, Rp2, R0, R1, R2)
+        #RESETS CONSTANTS HERE, the GLOBAL ROTATION MATRIX
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'tilt=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # tilt angle
+    if test:
+        tilt = float(test.group(1))
+        if not (tilt > -180) or not (tilt < 180):
+            errMsg(
+                "TILT=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
+            ) % tilt
             sys.exit(1)
-        if arg == "-":
-            break  #fixme
-        test = re.match(r'pan=(\d+\.\d+)', arg, flags=re.IGNORECASE)
-        if test:
-            pan = float(test.group(1))
-            if (pan > 180) or (pan < -180):
-                errMsg("PAN=%s MUST BE GREATER THAN -180 AND LESS THAN +180"
-                ) % pan
-                sys.exit(1)
-            panang = numpy.multiply(numpy.pi, numpy.divide(pan, 180))
-            sinpan = numpy.sin(panang)
-            cospan = numpy.cos(panang)
-            Rp0 = [cospan, 0, sinpan]
-            Rp1 = [0, 1, 0]
-            Rp2 = [-sinpan, 0, cospan]
-            # do matrix multiply to get new rotation matrix
-            newmat = matmul3(Rp0, Rp1, Rp2, R0, R1, R2)
-            #RESETS CONSTANTS HERE, the GLOBAL ROTATION MATRIX
-            R0 = newmat[0]
-            R1 = newmat[1]
-            R2 = newmat[2]
-        test = re.match(r'tilt=(\d+\.\d+)', arg, flags=re.IGNORECASE) # tilt angle
-        if test:
-            tilt = float(test.group(1))
-            if not (tilt > -180) or not (tilt < 180):
-                errMsg("TILT=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
-                ) % tilt
-                sys.exit(1)
-            tiltang = numpy.multiply(numpy.py, numpy.divide(tilt, 180))
-            sintilt = numpy.sin(tiltang)
-            costilt = numpy.cos(tiltang)
-            Rt0 = [1, 0, 0]
-            Rt1 = [0, costilt, sintilt]
-            Rt2 = [0, sintiltm, costilt]
-            # do matrix multiply to get new rotation matrix
-            newmat = matmul3(Rt0, Rt1, Rt2, R0, R1, R2)
-            #again resets the matrix ... hmm
-            R0 = newmat[0]
-            R1 = newmat[1]
-            R2 = newmat[2]
-        test = re.match(r'roll=(\d+\.\d+)', arg, flags=re.IGNORECASE) # roll angle
-        if test:
-            roll = float(test.group(1))
-            if not (roll > -180) or not (roll < 180):
-                errMsg("ROLL=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
-                ) % roll
+        tiltang = numpy.multiply(numpy.py, numpy.divide(tilt, 180))
+        sintilt = numpy.sin(tiltang)
+        costilt = numpy.cos(tiltang)
+        Rt0 = [1, 0, 0]
+        Rt1 = [0, costilt, sintilt]
+        Rt2 = [0, -sintilt, costilt]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rt0, Rt1, Rt2, R0, R1, R2)
+        #again resets the matrix ... hmm
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'roll=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # roll angle
+    if test:
+        roll = float(test.group(1))
+        if not (roll > -180) or not (roll < 180):
+            errMsg(
+                "ROLL=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
+            ) % roll
 
-            rollang = numpy.multiply(numpy.py, numpy.divide(roll, 180))
-            sinroll = numpy.sin(rollang)
-            cosroll = numpy.cos(rollang)
-            Rr0 = [cosroll, sinroll, 0]
-            Rr1 = [-sinroll, cosroll, 0]
-            Rr2 = [0, 0, 1]
-            # do matrix multiply to get new rotation matrix
-            newmat = matmul3(Rr0, Rr1, Rr2, R0, R1, R2)
-            R0 = newmat[0]
-            R1 = newmat[1]
-            R2 = newmat[2]
-        test = re.match(r'pef=(\d+\.\d+)', arg, flags=re.IGNORECASE) # pef angle
-        if test:
-            pef = float(test.group(1))
-            if not (pef > 0) or not (pef < numpy.pi):
-                errMsg("PEF=%s MUST BE LESS THAN PI AND GREATER THAN 0"
-                ) % pef
-        test = re.match(r'idx=(\d+\.\d+)', arg, flags=re.IGNORECASE) # idx angle
-        if test:
-            idx = float(test.group(1))
-            if not (idx >= 0):
-                errMsg("IDX=%s MUST BE 0 or GREATER"
-                ) % idx
-            di = idx
-        test = re.match(r'idy=(\d+\.\d+)', arg, flags=re.IGNORECASE)
-        if test:
-            idy = float(test.group(1))
-            if not (idy >= 0):
-                errMsg("IDY=%s MUST BE 0 or GREATER"
-                ) % idy
-            dj = idy
-        test = re.match(r'odx=(\d+\.\d+)', arg, flags=re.IGNORECASE)
-        if test:
-            odx = float(test.group(1))
-            if not (odx >= 0):
-                errMsg("ODX=%s MUST BE 0 or GREATER"
-                ) % odx
-            du = odx
-        test = re.match(r'ody=(\d+\.\d+)', arg, flags=re.IGNORECASE)
-        if test:
-            ody = float(test.group(1))
-            if not (ody > 0):
-                errMsg("ODY=%s MUST BE 0 or GREATER"
-                ) % ody
-            dv = ody
-        test = re.match(r'zoom=(\d+\.\d+)', arg, flags=re.IGNORECASE)
-        if test:
-            zoom = float(test.group(1))
-            if not (zoom < -1 or zoom > 1):
-                errMsg("ODY=%s MUST BE GREATER THAN 1 OR LESS THAN -1"
-                ) % zoom
-        test = re.match(r'bgcolor=([^ ]+)', arg, flags=re.IGNORECASE)
-        if test:
-            bgcolor = test.group(1)
-        test = re.match(r'skycolor=([^ ]+)', arg, flags=re.IGNORECASE)
-        if test:
-            skycolor = test.group(1)
-        test = re.match(r'vp=([^ ]+)', arg, flags=re.IGNORECASE)
-        if test:
-            vp = test.group(1).lower()
-            if vp not in [
-                    "background", "dither", "edge",
-                    "mirror", "random", "tile", "transparent"
-            ]:
-                errMsg("VP=%s IS NOT A VALID VALUE" % (vp))
-        test = re.match(r'auto=([^ ]+)', arg, flags=re.IGNORECASE)
-        if test:
-            auto = test.group(1).lower()
-            if auto not in [
-                    "c", "zc", "out"
-            ]:
-                errMsg("AUTO=%s IS NOT A VALID VALUE" % (auto))
-    #
-	# get infile and outfile
-    infile = sys.argv[-2]
-    outfile = sys.argv[-1]
-    for a in [infile, outfile]:
-        if re.match(r'.*=.*', a):
-            errMsg("%s is not a valid filepath" % a)
+        rollang = numpy.multiply(numpy.py, numpy.divide(roll, 180))
+        sinroll = numpy.sin(rollang)
+        cosroll = numpy.cos(rollang)
+        Rr0 = [cosroll, sinroll, 0]
+        Rr1 = [-sinroll, cosroll, 0]
+        Rr2 = [0, 0, 1]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rr0, Rr1, Rr2, R0, R1, R2)
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'pef=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # pef angle
+    if test:
+        pef = float(test.group(1))
+        if not (pef > 0) or not (pef < numpy.pi):
+            errMsg(
+                "PEF=%s MUST BE LESS THAN PI AND GREATER THAN 0"
+            ) % pef
+    test = re.match(r'idx=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # idx angle
+    if test:
+        idx = float(test.group(1))
+        if not (idx >= 0):
+            errMsg(
+                "IDX=%s MUST BE 0 or GREATER"
+            ) % idx
+        di = idx
+    test = re.match(r'idy=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        idy = float(test.group(1))
+        if not (idy >= 0):
+            errMsg(
+                "IDY=%s MUST BE 0 or GREATER"
+            ) % idy
+        dj = idy
+    test = re.match(r'odx=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        odx = float(test.group(1))
+        if not (odx >= 0):
+            errMsg(
+                "ODX=%s MUST BE 0 or GREATER"
+            ) % odx
+        du = odx
+    test = re.match(r'ody=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        ody = float(test.group(1))
+        if not (ody > 0):
+            errMsg(
+                "ODY=%s MUST BE 0 or GREATER"
+            ) % ody
+        dv = ody
+    test = re.match(r'zoom=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        zoom = float(test.group(1))
+        if not (zoom < -1 or zoom > 1):
+            errMsg(
+                "ODY=%s MUST BE GREATER THAN 1 OR LESS THAN -1"
+            ) % zoom
+    test = re.match(r'bgcolor=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        bgcolor = test.group(1)
+    test = re.match(r'skycolor=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        skycolor = test.group(1)
+    test = re.match(r'vp=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        vp = test.group(1).lower()
+        if vp not in [
+                "background", "dither", "edge",
+                "mirror", "random", "tile", "transparent"
+        ]:
+            errMsg("VP=%s IS NOT A VALID VALUE" % (vp))
+    test = re.match(r'auto=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        auto = test.group(1).lower()
+        if auto not in [
+                "c", "zc", "out"
+        ]:
+            errMsg("AUTO=%s IS NOT A VALID VALUE" % (auto))
+#
+# get infile and outfile
+infile = sys.argv[-2]
+outfile = sys.argv[-1]
+for a in [infile, outfile]:
+    if re.match(r'.*=.*', a):
+        errMsg("%s is not a valid filepath" % a)
 
 
-    # setup temporary images and auto delete upon exit
-    # use mpc/cache to hold input image temporarily in memory
-    pid = os.getpid()
+# setup temporary images and auto delete upon exit
+# use mpc/cache to hold input image temporarily in memory
+pid = os.getpid()
 
-    tmpA = "%s/3Drotate_%s.mpc" % (directory, pid)
-    tmpB = "%s/3Drotate_%s.cache" % (directory, pid)
-    os.unlink(tmpA,tmpB)
-    os.unlink(tmpA,tmpB)
+tmpA = "%s/3Drotate_%s.mpc" % (directory, pid)
+tmpB = "%s/3Drotate_%s.cache" % (directory, pid)
+os.unlink(tmpA, tmpB)
+os.unlink(tmpA, tmpB)
 
-    # test that infile provided
-    if not os.path.exists(infile):
-        errMsg("NO INPUT FILE SPECIFIED")
-    # test that outfile provided
-    if not outfile:
-        errMsg("NO OUTPUT FILE SPECIFIED")
+# test that infile provided
+if not os.path.exists(infile):
+    errMsg("NO INPUT FILE SPECIFIED")
+# test that outfile provided
+if not outfile:
+    errMsg("NO OUTPUT FILE SPECIFIED")
 
-    cmd = ['convert', '-quiet', '-regard-warnings', infile, '+repage', tmpA]
-    call_cmd(cmd)
-    if not pef:
-        pef = 1
+cmd = ['convert', '-quiet', '-regard-warnings', infile, '+repage', tmpA]
+call_cmd(cmd)
+if not pef:
+    pef = 1
 
-    # get input image width and height
-    imagesize()
-    maxwidth = width - 1
-    maxheight = height - 1
+# get input image width and height
+width, height = imagesize()
+maxwidth = width - 1
+maxheight = height - 1
 
-    # deal with auto adjustments to values
-    if auto == "zc":
-        du = 0
-        dv = 0
-        zoom = 1
-    elif auto == "c":
-        du = 0
-        dv = 0
+# deal with auto adjustments to values
+if auto == "zc":
+    du = 0
+    dv = 0
+    zoom = 1
+elif auto == "c":
+    du = 0
+    dv = 0
 
-    # convert offsets of rotation point to relative to pixel 0,0
-    di = numpy.divide(numpy.add(di, numpy.divide(numpy.subtract(width,1), 2)),1)
-    dj = numpy.divide(
-            numpy.add(dj, numpy.divide(numpy.subtract(height, 1), 2)), 1)
-    du = numpy.divide(du, 1)
-    dv = numpy.divide(dv, 1)
+# convert offsets of rotation point to relative to pixel 0,0
+di = numpy.divide(numpy.add(di, numpy.divide(numpy.subtract(width, 1), 2)), 1)
+dj = numpy.divide(
+        numpy.add(dj, numpy.divide(numpy.subtract(height, 1), 2)), 1)
+du = numpy.divide(du, 1)
+dv = numpy.divide(dv, 1)
 
-    # convert zoom to scale factors
-    if zoom >= 1:
-        sx = numpy.divide(1, zoom)
-        sy = sx
-    elif zoom <= -1:
-        sx = numpy.divide(-zoom, 1)
-        sy = sx
+# convert zoom to scale factors
+if zoom >= 1:
+    sx = numpy.divide(1, zoom)
+    sy = sx
+elif zoom <= -1:
+    sx = numpy.divide(-zoom, 1)
+    sy = sx
 
-    #{{{explanation
-    # Consider the picture placed on the Z=0 plane and the camera a distance
-    # Zc=f above the picture plane looking straight down at the image center.
-    # Now the perspective equations (in 3-D) are defined as (x,y,f) = M (X',Y',Z'),
-    # where the camera orientation matrix M is the identity matrix but with M22=-1
-    # because the camera is looking straight down along -Z.
-    # Thus a reflection transformation relative to the ground plane coordinates.
-    # Let the camera position Zc=f=(sqrt(ins*ins + inl*inl)) / ( 2 tan(fov/2) )
-    # Now we want to rotate the ground points corresponding to the picture corners.
-    # The basic rotation is (X',Y',Z') = R (X,Y,0), where R is the rotation matrix
-    # involving pan, tilt and roll.
-    # But we need to convert (X,Y,0) to (X,Y,1) and also to offset for Zc=f
-    # First we note that (X,Y,0) = (X,Y,1) - (0,0,1)
-    # Thus the equation becomes (x,y,f) = M {R [(X,Y,1) - (0,0,1)] - (0,0,Zc)} = MT (X,Y,1)
-    # But R [(X,Y,1) - (0,0,1)] = R [II (X,Y,1) - S (X,Y,1)] = R (II-S) (X,Y,1), where
-    # II is the identity matrix and S is an all zero matrix except for S22=1.
-    # Thus (II-S) is the identity matrix with I22=0 and
-    # RR = R (II-S) is just R with the third column all zeros.
-    # Thus we get (x,y,f) = M {RR (X,Y,1) - (0,0,Zc)}.
-    # But M {RR (X,Y,1) - (0,0,Zc)} = M {RR(X,Y,1) - D (X,Y,1)}, where
-    # D is an all zero matrix with D22 = Zc = f.
-    # So that we get M (RR-D) (X,Y,1) = MT (X,Y,1), where
-    # where T is just R with the third column (0,0,-f), i.e. T02=0, T12=0, T22=-f
-    # But we need to allow for scaling and offset of the output coordinates and
-    # conversion from (x,y,f) to (u,v,1)=O and conversion of input coordinates
-    # from (X,Y,1) to (i,j,1)=I.
-    # Thus the forward transformation becomes AO=MTBI or O=A'MTBI or O=PI,
-    # where prime means inverse.
-    # However, to do the scaling of the output correctly, need to offset by the input
-    # plus output offsets, then scale, which is all put into A'.
-    # Thus the forward transformation becomes AO=MTBI or O=A'MTBI where A'=Ai
-    # but we will merge A'M into Aim
-    # Thus the inverse transform becomes
-    # I=QO where Q=P'
-    # A=output scaling, offset and conversion matrix
-    # B=input offset and conversion matrix (scaling only needs to be done in one place)
-    # M=camera orientation matrix
-    # R=image rotation matrix Rroll Rtilt Rpan
-    # T=matrix that is R but R33 offset by f + 1
-    # O=output coords vector (i,j,1)
-    # I=input coords vector (u,v,1)=(is,il,1)
-    # P=forward perspective transformation matrix
-    # Q=inverse perspective transformation matrix
-    #
-    # For a 35 mm camera whose film format is 36mm wide and 24mm tall, when the focal length
-    # is equal to the diagonal, the field of view is 53.13 degrees and this is
-    # considered a normal view equivalent to the human eye.
-    # See http://www.panoramafactory.com/equiv35/equiv35.html
-    # Max limit on dfov is 180 degrees (pef=3.19) where get single line like looking at picture on edge.
-    # Above this limit the picture becomes like the angles get reversed.
-    # Min limit on dfov seems to be slightly greater than zero degrees.
-    # Practical limits on dfov depend upon orientation angles.
-    # For tilt=45, this is about 2.5 dfov (pef=2.5). Above this, some parts of the picture
-    # that are cut off at the bottom, get wrapped and stretched in the 'sky'.
-    #}}}
+#{{{explanation
+# Consider the picture placed on the Z=0 plane and the camera a distance
+# Zc=f above the picture plane looking straight down at the image center.
+# Now the perspective equations (in 3-D) are defined as (x,y,f) = M (X',Y',Z'),
+# where the camera orientation matrix M is the identity matrix but with M22=-1
+# because the camera is looking straight down along -Z.
+# Thus a reflection transformation relative to the ground plane coordinates.
+# Let the camera position Zc=f=(sqrt(ins*ins + inl*inl)) / ( 2 tan(fov/2) )
+# Now we want to rotate the ground points corresponding to the picture corners.
+# The basic rotation is (X',Y',Z') = R (X,Y,0), where R is the rotation matrix
+# involving pan, tilt and roll.
+# But we need to convert (X,Y,0) to (X,Y,1) and also to offset for Zc=f
+# First we note that (X,Y,0) = (X,Y,1) - (0,0,1)
+# Thus the equation becomes
+# (x,y,f) = M {R [(X,Y,1) - (0,0,1)] - (0,0,Zc)} = MT (X,Y,1)
+# But R [(X,Y,1) - (0,0,1)] = R [II (X,Y,1) - S (X,Y,1)] = R (II-S) (X,Y,1),
+# where II is the identity matrix and S is an all zero matrix except for S22=1.
+# Thus (II-S) is the identity matrix with I22=0 and
+# RR = R (II-S) is just R with the third column all zeros.
+# Thus we get (x,y,f) = M {RR (X,Y,1) - (0,0,Zc)}.
+# But M {RR (X,Y,1) - (0,0,Zc)} = M {RR(X,Y,1) - D (X,Y,1)}, where
+# D is an all zero matrix with D22 = Zc = f.
+# So that we get M (RR-D) (X,Y,1) = MT (X,Y,1), where
+# where T is just R with the third column (0,0,-f), i.e. T02=0, T12=0, T22=-f
+# But we need to allow for scaling and offset of the output coordinates and
+# conversion from (x,y,f) to (u,v,1)=O and conversion of input coordinates
+# from (X,Y,1) to (i,j,1)=I.
+# Thus the forward transformation becomes AO=MTBI or O=A'MTBI or O=PI,
+# where prime means inverse.
+# However, to do the scaling of the output correctly, need to offset by the
+# input plus output offsets, then scale, which is all put into A'.
+# Thus the forward transformation becomes AO=MTBI or O=A'MTBI where A'=Ai
+# but we will merge A'M into Aim
+# Thus the inverse transform becomes
+# I=QO where Q=P'
+# A=output scaling, offset and conversion matrix
+# B=input offset and conversion matrix
+# (scaling only needs to be done in one place)
+# M=camera orientation matrix
+# R=image rotation matrix Rroll Rtilt Rpan
+# T=matrix that is R but R33 offset by f + 1
+# O=output coords vector (i,j,1)
+# I=input coords vector (u,v,1)=(is,il,1)
+# P=forward perspective transformation matrix
+# Q=inverse perspective transformation matrix
+#
+# For a 35 mm camera whose film format is 36mm wide
+# and 24mm tall, when the focal length
+# is equal to the diagonal, the field of view is 53.13 degrees and this is
+# considered a normal view equivalent to the human eye.
+# See http://www.panoramafactory.com/equiv35/equiv35.html
+# Max limit on dfov is 180 degrees (pef=3.19)
+# where get single line like looking at picture on edge.
+# Above this limit the picture becomes like the angles get reversed.
+# Min limit on dfov seems to be slightly greater than zero degrees.
+# Practical limits on dfov depend upon orientation angles.
+# For tilt=45, this is about 2.5 dfov (pef=2.5).
+# Above this, some parts of the picture
+# that are cut off at the bottom, get wrapped and stretched in the 'sky'.
+#}}}
 
-    dfov = numpy.divide(
-            numpy.multpily(180, numpy.arctan(36/24)), numpy.pi)
-    pfact = 1
-    if pef == 0: #fixme
-        pfact = numpy.divide(0.01, dfov)
-    else:
-        pfact = pef
+dfov = numpy.divide(
+        numpy.multpily(180, numpy.arctan(36/24)), numpy.pi)
+pfact = 1
+if pef == 0:  # FIXME
+    pfact = numpy.divide(0.01, dfov)
+else:
+    pfact = pef
 
-    #compute new field of view based upon pef (pfact)
-    dfov = numpy.multiply(pfact, dfov)
-    dfov2 = numpy.divide(dfov, 2)
-    arg = numpy.multiply(numpy.pi, numpy.divide(dfov, 180))
-    sfov = numpy.sin(arg)
-    cfov = numpy.cos(arg)
-    tfov = numpy.divide(sfov, cfov)
+#compute new field of view based upon pef (pfact)
+dfov = numpy.multiply(pfact, dfov)
+dfov2 = numpy.divide(dfov, 2)
+arg = numpy.multiply(numpy.pi, numpy.divide(dfov, 180))
+sfov = numpy.sin(arg)
+cfov = numpy.cos(arg)
+tfov = numpy.divide(sfov, cfov)
 
 
-    # calculate focal length in same units as wall (picture) using dfov
-    diag = numpy.sqrt(
-            numpy.multiply(width, width), numpy.multiply(height, height))
-    focal = numpy.divide(diag, numpy.multiply(2, tfov))
+# calculate focal length in same units as wall (picture) using dfov
+diag = numpy.sqrt(
+        numpy.multiply(width, width), numpy.multiply(height, height))
+focal = numpy.divide(diag, numpy.multiply(2, tfov))
 
-    # calculate forward transform matrix Q
+# calculate forward transform matrix Q
 
-    # define the input offset and conversion matrix
-    dim = -di
-    B0 = [1, 0, dim]
-    B1 = [0, -1, dj]
-    B2 = [0, 0, 1]
+# define the input offset and conversion matrix
+dim = -di
+B0 = [1, 0, dim]
+B1 = [0, -1, dj]
+B2 = [0, 0, 1]
 
-    # define the output scaling, offset and conversion matrix inverse Ai and merge with M
-    # to become Aim
-    #A0=($sx 0 $sx*(-$du-$di))
-    #A1=(0 -$sy $sy*($dv+$dj))
-    #A2=(0 0 -$focal)
-    #M0=(1 0 0)
-    #M1=(0 1 0)
-    #M2=(0 0 -1)
-    aim00 = numpy.divide(1, sx)
-    aim02 = -numpy.divide(
-            numpy.multiply(sx, numpy.add(di, du)), numpy.multiply(sx, focal))
-    aim11 = numpy.divide(-1, sy)
-    aim12 = -numpy.divide(
-            numpy.multiply(sy, numpy.add(dj, dv)), numpy.multiply(sy, focal))
-    aim22 = -numpy.divide(1, focal)
-    Aim0 = [aim00, 0, aim02]
-    Aim1 = [0, aim11, aim12]
-    Aim2 = [0, 0, aim22]
+# define the output scaling, offset and conversion
+# matrix inverse Ai and merge with M
+# to become Aim
+#A0=($sx 0 $sx*(-$du-$di))
+#A1=(0 -$sy $sy*($dv+$dj))
+#A2=(0 0 -$focal)
+#M0=(1 0 0)
+#M1=(0 1 0)
+#M2=(0 0 -1)
+aim00 = numpy.divide(1, sx)
+aim02 = -numpy.divide(
+        numpy.multiply(sx, numpy.add(di, du)), numpy.multiply(sx, focal))
+aim11 = numpy.divide(-1, sy)
+aim12 = -numpy.divide(
+        numpy.multiply(sy, numpy.add(dj, dv)), numpy.multiply(sy, focal))
+aim22 = -numpy.divide(1, focal)
+Aim0 = [aim00, 0, aim02]
+Aim1 = [0, aim11, aim12]
+Aim2 = [0, 0, aim22]
 
-    # now do successive matrix multiplies
-    # from right towards left of main equation P=A'RB
+# now do successive matrix multiplies
+# from right towards left of main equation P=A'RB
 
-    # convert R to T by setting T02=T12=0 and T22=-f
-    T0 = [R0[0], R0[1], 0]
-    T1 = [R1[0], R1[1], 0]
-    T2 = [R2[0], R2[1], -focal]
+# convert R to T by setting T02=T12=0 and T22=-f
+T0 = [R0[0], R0[1], 0]
+T1 = [R1[0], R1[1], 0]
+T2 = [R2[0], R2[1], -focal]
 
-    # multiply T x B = P
-    p_mat = matmul3(T0, T1, T2, B0, B1, B2)
+# multiply T x B = P
+p_mat = matmul3(T0, T1, T2, B0, B1, B2)
 
-    # multiply Aim x P = P
-    newmat = matmul3(Aim0, Aim1, Aim2, p_mat[0], p_mat[1], p_mat[2])
+# multiply Aim x P = P
+newmat = matmul3(Aim0, Aim1, Aim2, p_mat[0], p_mat[1], p_mat[2])
 
-    #    # the resulting P matrix is now
-    #    the perspective coefficients for the inverse transformation
-    #    P00= newmat[0][0]
-    #    P01= newmat[0][1]
-    #    P02= newmat[0][2]
-    #    P10= newmat[1][0]
-    #    P11= newmat[1][1]
-    #    P12= newmat[1][2]
-    #    P20= newmat[2][0]
-    #    P21= newmat[2][1]
-    #    P22= newmat[2][2]
+#    # the resulting P matrix is now
+#    the perspective coefficients for the inverse transformation
+#    P00= newmat[0][0]
+#    P01= newmat[0][1]
+#    P02= newmat[0][2]
+#    P10= newmat[1][0]
+#    P11= newmat[1][1]
+#    P12= newmat[1][2]
+#    P20= newmat[2][0]
+#    P21= newmat[2][1]
+#    P22= newmat[2][2]
 
-    # project input corners to output domain
-    #echo "UL"
-    i = 0
-    j = 0
-    u1, v1 = forwardProject(newmat, i, j) #using Aim matrix and p_mat
+# project input corners to output domain
+#echo "UL"
+i = 0
+j = 0
+u1, v1 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
 
-    i = maxwidth
-    j = 0
-    u2, v2 = forwardProject(newmat, i, j) #using Aim matrix and p_mat
+i = maxwidth
+j = 0
+u2, v2 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
 
-    i = maxwidth
-    j = maxheight
+i = maxwidth
+j = maxheight
 
-    u3, v3 = forwardProject(newmat, i, j) #using Aim matrix and p_mat
-    i = 0
-    j = maxheight
+u3, v3 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
+i = 0
+j = maxheight
 
-    u4, v4 = forwardProject(newmat, i, j) #using Aim matrix and p_mat
+u4, v4 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
 
-    # deal with adjustments for auto settings
-    # first get the bounding box dimensions
-    uArr = [u1, u2, u3, u4]
-    vArr = [v1, v2, v3, v4]
-    umin = 1000000
-    umax = -1000000
-    vmin = 1000000
-    vmax = -1000000
-    for index in range(0, 4):
-        if uArr[index] < umin:
-            umin = uArr[index]
-        if uArr[index] > umax:
-            umax = uArr[index]
-        if vArr[index] < vmin:
-            vmin=vArr[index]
-        if vArr[index] > vmax:
-            vmax=vArr[index]
-    delu = numpy.add(numpy.subtract(umax, umin), 1)
-    delv = numpy.add(numpy.subtract(vmax, vmin), 1)
-    if auto == "c":
-        offsetu = `echo "scale = 10; ($width - $delu) / 2" | bc`
-        offsetv = `echo "scale = 10; ($height - $delv) / 2" | bc`
-        u1 = `echo "scale = 0; $offsetu + ($u1 - $umin)" | bc`
-        v1 = `echo "scale = 0; $offsetv + ($v1 - $vmin)" | bc`
-        u2 = `echo "scale = 0; $offsetu + ($u2 - $umin)" | bc`
-        v2 = `echo "scale = 0; $offsetv + ($v2 - $vmin)" | bc`
-        u3 = `echo "scale = 0; $offsetu + ($u3 - $umin)" | bc`
-        v3 = `echo "scale = 0; $offsetv + ($v3 - $vmin)" | bc`
-        u4 = `echo "scale = 0; $offsetu + ($u4 - $umin)" | bc`
-        v4 = `echo "scale = 0; $offsetv + ($v4 - $vmin)" | bc`
-elif [ "$auto" = "zc" ]
-	then
-	if [ `echo "$delu > $delv" | bc` -eq 1 ]
-		then
-		del=$delu
-		offsetu=0
-		offsetv=`echo "scale=10; ($height - ($delv * $width / $delu)) / 2" | bc`
-	else
-		del=$delv
-		offsetu=`echo "scale=10; ($width - ($delu * $height / $delv)) / 2" | bc`
-		offsetv=0
-	fi
-	u1 = numpy.add(
-        offsetu,
-        numpy.multpily(numpy.subtract(u1, umin), numpy.divide(width, del)))
-	v1 = numpy.add(
-        offsetv,
-        numpy.multpily(numpy.subtract(v1, vmin), numpy.divide(height, del)))
-	u2 = numpy.add(
-        offsetu,
-        numpy.multpily(numpy.subtract(u2, umin), numpy.divide(width, del)))
-	v2 = numpy.add(
-        offsetv,
-        numpy.multpily(numpy.subtract(v2, vmin), numpy.divide(height, del)))
-	u3 = numpy.add(
-        offsetu,
-        numpy.multpily(numpy.subtract(u3, umin), numpy.divide(width, del)))
-	v3 = numpy.add(
-        offsetv,
-        numpy.multpily(numpy.subtract(v3, vmin), numpy.divide(height, del)))
-	u4 = numpy.add(
-        offsetu,
-        numpy.multpily(numpy.subtract(u4, umin), numpy.divide(width, del)))
-	v4 = numpy.add(
-        offsetv,
-        numpy.multpily(numpy.subtract(v4, vmin), numpy.divide(height, del)))
-fi
+# deal with adjustments for auto settings
+# first get the bounding box dimensions
+uArr = [u1, u2, u3, u4]
+vArr = [v1, v2, v3, v4]
+umin = 1000000
+umax = -1000000
+vmin = 1000000
+vmax = -1000000
+for index in range(0, 4):
+    if uArr[index] < umin:
+        umin = uArr[index]
+    if uArr[index] > umax:
+        umax = uArr[index]
+    if vArr[index] < vmin:
+        vmin = vArr[index]
+    if vArr[index] > vmax:
+        vmax = vArr[index]
+delu = numpy.add(numpy.subtract(umax, umin), 1)
+delv = numpy.add(numpy.subtract(vmax, vmin), 1)
+if auto == "c":
+    offsetu = numpy.divide(numpy.subtract(width, delu), 2)
+    offsetv = numpy.divide(numpy.subtract(height, delv), 2)
+    #FIXME may need to cast as int variables below (u1-v4)
+    u1 = numpy.add(offsetu, numpy.subtract(u1, umin))
+    v1 = numpy.add(offsetv, numpy.subtract(v1, vmin))
+    u2 = numpy.add(offsetu, numpy.subtract(u2, umin))
+    v2 = numpy.add(offsetv, numpy.subtract(v2, vmin))
+    u3 = numpy.add(offsetu, numpy.subtract(u3, umin))
+    v3 = numpy.add(offsetv, numpy.subtract(v3, vmin))
+    u4 = numpy.add(offsetu, numpy.subtract(u4, umin))
+    v4 = numpy.add(offsetv, numpy.subtract(v4, vmin))
+elif auto == "zc":
+    if delu > delv:
+        _del = delu
+        offsetu = 0
+        offsetv = numpy.divide(
+            numpy.subtract(
+                height, numpy.divide(numpy.multpily(delv, width), delu)), 2)
+    else:
+        _del = delv
+        offsetu = numpy.divide(
+            numpy.subtract(
+                width, numpy.divide(numpy.multpily(delu, height), delv)), 2)
+        offsetv = 0
+u1 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u1, umin), numpy.divide(width, _del)))
+v1 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v1, vmin), numpy.divide(height, _del)))
+u2 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u2, umin), numpy.divide(width, _del)))
+v2 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v2, vmin), numpy.divide(height, _del)))
+u3 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u3, umin), numpy.divide(width, _del)))
+v3 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v3, vmin), numpy.divide(height, _del)))
+u4 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u4, umin), numpy.divide(width, _del)))
+v4 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v4, vmin), numpy.divide(height, _del)))
 #
 # now do the perspective distort
-if [ "$auto" = "out" ]
-	then
-	distort="+distort"
-else
-	distort="-distort"
-fi
+if auto == "out":
+    distort = "+distort"
+else:
+    distort = "-distort"
 
-im_version=`convert -list configure | \
-	sed '/^LIB_VERSION_NUMBER /!d; s//,/;  s/,/,0/g;  s/,0*\([0-9][0-9]\)/\1/g' | head -n 1`
-if [ "$im_version" -lt "06030600" ]
-	then
-	convert $tmpA -virtual-pixel $vp -background $bgcolor \
-	-mattecolor $skycolor $distort Perspective \
-	"0,0 $maxwidth,0 $maxwidth,$maxheight 0,$maxheight  $u1,$v1 $u2,$v2 $u3,$v3 $u4,$v4" $outfile
-else
-	convert $tmpA -virtual-pixel $vp -background $bgcolor \
-	-mattecolor $skycolor $distort Perspective \
-	"0,0 $u1,$v1   $maxwidth,0 $u2,$v2   $maxwidth,$maxheight $u3,$v3   0,$maxheight $u4,$v4" $outfile
-fi
-exit 0
+#im_version=`convert -list configure | \
+#	sed '/^LIB_VERSION_NUMBER /!d; s//,/;  s/,/,0/g;  s/,0*\([0-9][0-9]\)/\1/g'\
+#	| head -n 1`
+#if [ "$im_version" -lt "06030600" ]
+#	then
+#	convert $tmpA -virtual-pixel $vp -background $bgcolor \
+#	-mattecolor $skycolor $distort Perspective \
+#	"0,0 $maxwidth,0 $maxwidth,$maxheight 0,$maxheight  \
+#	$u1,$v1 $u2,$v2 $u3,$v3 $u4,$v4" $outfile
+#else
+cmd = [
+    "convert", tmpA, "-virtual-pixel", vp, "-background", bgcolor,
+    "-mattecolor", skycolor, distort, "Perspective",
+    "0,0", "%s,%s" % (u1, v1),
+    "%s,0" % (maxwidth), "%s,%s" % (u2, v2),
+    "%s,%s" % (maxwidth, maxheight),
+    "%s,%s" % (u3, v3),
+    "%s,%s" % (0, maxheight),
+    "%s,%s" % (u4, v4),
+    outfile
+]
+call_cmd(cmd)
diff --git a/bin/3Drotate.py.bak b/bin/3Drotate.py.bak
new file mode 100755
index 0000000..a94a7a9
--- /dev/null
+++ b/bin/3Drotate.py.bak
@@ -0,0 +1,716 @@
+#!/usr/bin/python2.7
+import re
+import numpy
+import sys
+import os
+from subprocess import Popen, PIPE
+
+#{{{ usage
+USAGE = """
+USAGE: 3Drotate option=value infile outfile
+USAGE: 3Drotate [-h or -help]
+
+OPTIONS:  any one or more
+
+pan      value     rotation about image vertical centerline;
+                   -180 to +180 (deg); default=0
+tilt     value     rotation about image horizontal centerline;
+                   -180 to +180 (deg); default=0
+roll     value     rotation about the image center;
+                   -180 to +180 (deg); default=0
+pef      value     perspective exaggeration factor;
+                   0 to 3.19; default=1
+idx  	   value     +/- pixel displacement in rotation point right/left
+                   in input from center; default=0
+idy      value     +/- pixel displacement in rotation point down/up
+                   in input from center; default=0
+odx      value     +/- pixel displacement in rotation point right/left
+                   in output from center; default=0
+ody      value     +/- pixel displacement in rotation point down/up
+                   in output from center; default=0
+zoom     value     output zoom factor; where value > 1 means zoom in
+                   and < -1 means zoom out; value=1 means no change
+bgcolor  value     the background color value; any valid IM image
+                   color specification (see -fill); default is black
+skycolor value     the sky color value; any valid IM image
+                   color specification (see -fill); default is black
+auto     c         center bounding box in output
+                   (odx and ody ignored)
+auto     zc        zoom to fill and center bounding box in output
+                   (odx, ody and zoom ignored)
+auto     out       creates an output image of size needed to hold
+                   the transformed image; (odx, ody and zoom ignored)
+vp       value     virtual-pixel method; any valid IM virtual-pixel method;
+                   default=background
+
+#
+
+NAME: 3DROTATE
+
+PURPOSE: To apply a perspective distortion to an image by
+providing rotation angles, zoom, offsets, background color,
+perspective exaggeration and auto zoom/centering.
+
+DESCRIPTION: 3DROTATE applies a perspective distortion to an image
+by providing any combination of three optional rotation angle:
+pan, tilt and roll with optional offsets and zoom and with an optional
+control of the perspective exaggeration. The image is treated as if it
+were painted on the Z=0 ground plane. The picture plane is then rotated
+and then perspectively projected to a camera located a distance equal to
+the focal length above the ground plane looking straight down along
+the -Z direction.
+
+
+ARGUMENTS:
+
+PAN is a rotation of the image about its vertical
+centerline -180 to +180 degrees. Positive rotations turn the
+right side of the image away from the viewer and the left side
+towards the viewer. Zero is no rotation. A PAN of +/- 180 deg
+achieves the same results as -flip.
+
+TILT is a rotation of the image about its horizontal
+centerline -180 to +180 degrees. Positive rotations turn the top
+of the image away from the viewer and the bottom towards the
+viewer. Zero is no rotation. A TILT of +/- 180 deg
+achieves the same results as -flop.
+
+ROLL (like image rotation) is a rotation in the plane of the
+the image -180 to +180 degrees. Positive values are clockwise
+and negative values are counter-clockwise. Zero is no rotation.
+A ROLL of any angle achieves the same results as -rotate.
+
+PAN, TILT and ROLL are order dependent. If all three are provided,
+then they will be done in whatever order specified.
+
+PEF is the perspective exaggeration factor. It ranges from 0 to 3.19.
+A normal perspective is achieved with the default of 1. As PEF is
+increased from 1, the perspective effect moves towards that of
+a wide angle lens (more distortion). If PEF is decreased from 1
+the perspective effect moves towards a telephoto lens (less
+distortion). PEF of 0.5 achieves an effect close to no perspective
+distortion. As pef gets gets larger than some value which depends
+upon the larger the pan, tilt and roll angles become, one reaches
+a point where some parts of the picture become so distorted that
+they wrap around and appear above the "horizon"
+
+IDX is the a pixel displacement of the rotation point in the input image
+from the image center. Positive values shift to the right along the
+sample direction; negative values shift to the left. The default=0
+corresponds to the image center.
+
+IDY is the a pixel displacement of the rotation point in the input image
+from the image center. Positive values shift to downward along the
+line direction; negative values shift upward. The default=0
+corresponds to the image center.
+
+ODX is the a pixel displacement from the center of the output image where
+one wants the corresponding input image rotation point to appear.
+Positive values shift to the right along the sample direction; negative
+values shift to the left. The default=0 corresponds to the output image center.
+
+ODY is the a pixel displacement from the center of the output image where
+one wants the corresponding input image rotation point to appear.
+Positive values shift downward along the sample direction; negative
+values shift upward. The default=0 corresponds to the output image center.
+
+ZOOM is the output image zoom factor. Values > 1 (zoomin) cause the image
+to appear closer; whereas values < 1 (zoomout) cause the image to
+appear further away.
+
+BGCOLOR is the color of the background to use to fill where the output image
+is outside the area of the perspective of the input image. See the IM function
+-fill for color specifications. Note that when using rgb(r,g,b), this must be
+enclosed in quotes after the equal sign.
+
+SKYCOLOR is the color to use in the 'sky' area above the perspective 'horizon'.
+See the IM function -fill for color specifications. Note that when using
+rgb(r,g,b), this must be enclosed in quotes after the equal sign.
+
+AUTO can be either c, zc or out. If auto is c, then the resulting perspective
+of the input image will have its bounding box centered in the output image
+whose size will be the same as the input image. If
+auto is zc, then the resulting perspective of the input image will have its
+bounding box zoomed to fill its largest dimension to match the size of the
+the input image and the other dimension will be centered in the output. If
+auto is out, then the output image will be made as large or as small as
+needed to just fill out the transformed input image. If any of these are
+present, then the arguments OSHIFTX, OSHIFTY are ignored.
+
+VP is the virtual-pixel method, which allows the image to be extended outside
+its bounds. For example, vp=background, then the background color is used to
+fill the area in the output image which is outside the perspective view of
+the input image. If vp=tile, then the perspective view will be tiled to fill
+the output image.
+
+NOTE: The output image size will be the same as the input image size due
+to current limitations on -distort Perspective.
+
+CAVEAT: No guarantee that this script will work on all platforms,
+nor that trapping of inconsistent parameters is complete and
+foolproof. Use At Your Own Risk.
+"""
+
+#}}}
+
+#{{{ defaults
+# set default value
+# rotation angles and rotation matrix
+pan = 0
+tilt = 0
+roll = 0
+R0 = [1, 0, 0]
+R1 = [0, 1, 0]
+R2 = [0, 0, 1]
+
+# scaling output only
+sx = 1
+sy = 1
+
+# offset du,dv = output; relative to center of image
+du = 0
+dv = 0
+
+# offset di,dj = input; relative to center of image
+di = 0
+dj = 0
+
+# perspective exaggeration factor
+pef = 1
+# zoom
+zoom = 1
+# background color
+bgcolor = "black"
+# sky color
+skycolor = "black"
+
+# virtual-pixel method
+vp = "background"
+
+# set directory for temporary files
+directory = "."    # suggestions are dir="." or dir="/tmp"
+
+pi = numpy.math.pi
+#}}}
+
+
+# function to do dot product of 2 three element vectors
+def dot(vec1, vec2):
+    return numpy.dot(vec1, vec2)
+
+
+# function to do 3x3 matrix multiply M x N where input
+# are rows of each matrix; M1 M2 M3 N1 N2 N3
+def matmul3(m0, m1, m2, n0, n1, n2):
+    return numpy.matmul(
+      [m0, m1, m2], [n0, n1, n2]
+    )
+
+
+# function to project points from input to output domain
+def forwardProject(mat, ii, jj):
+    #takes in a 3x3 matrix and scalars ii and jj
+    #returns two scalars, uu and vv
+    numu = mat[0][0] * ii + mat[0][1] * jj + mat[0][2]
+    numv = mat[1][0] * ii + mat[1][1] * jj + mat[1][2]
+    den = mat[2][0] * ii + mat[2][1] * jj + mat[2][2]
+    uu = numu / den
+    vv = numv / den
+    return [uu, vv]
+
+
+def errMsg(s):
+    sys.stderr.write("%s\n")
+    sys.exit(1)
+
+
+def usage():
+    sys.stderr.write(USAGE)
+    sys.exit(1)
+
+
+def inverseProject(mat, uu, vv):
+    #note this function is more exact in python version
+    numi = (mat[0][0] * uu) + (mat[0][1] * vv) + mat[0][2]
+    numj = (mat[1][0] * uu) + (mat[1][1] * vv) + mat[1][2]
+    den = (mat[2][0] * uu) + (mat[2][1] * vv) + mat[2][2]
+    ii = numi / den
+    jj = numj / den
+    #FIXME chop off decimal for it to be like bash version
+    return [ii, jj]
+
+
+# function to invert a 3 x 3 matrix using method of adjoint
+# inverse is the transpose of the matrix of cofactors divided by the determinant
+def M3inverse(mat):
+    return numpy.linalg.inv(mat)
+
+
+def call_cmd(s):
+    cmd = s.split(" ")
+    process = Popen(cmd, stdout=PIPE, stderr=PIPE)
+    out, err = process.communicate()
+    errcode = process.returncode
+    if errcode > 0:
+        errMsg(err)
+    return out
+
+
+# get input image size
+def imagesize(filepath):
+    width = int(call_cmd("identify -format %w %s" % (filepath)))
+    height = int(call_cmd("identify -format %h %s" % (filepath)))
+    return width, height
+
+
+# test for correct number of arguments and get values
+len_args = len(sys.argv)
+
+if len_args == 1:
+    usage()
+elif (len_args > 15):
+    errMsg("--- TOO MANY ARGUMENTS WERE PROVIDED ---")
+
+for arg in sys.argv:
+    if arg == "-h":
+        usage()
+    if arg == "-":
+        break  # FIXME
+    test = re.match(r'pan=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        pan = float(test.group(1))
+        if (pan > 180) or (pan < -180):
+            errMsg(
+                "PAN=%s MUST BE GREATER THAN -180 AND LESS THAN +180"
+            ) % pan
+            sys.exit(1)
+        panang = numpy.multiply(numpy.pi, numpy.divide(pan, 180))
+        sinpan = numpy.sin(panang)
+        cospan = numpy.cos(panang)
+        Rp0 = [cospan, 0, sinpan]
+        Rp1 = [0, 1, 0]
+        Rp2 = [-sinpan, 0, cospan]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rp0, Rp1, Rp2, R0, R1, R2)
+        #RESETS CONSTANTS HERE, the GLOBAL ROTATION MATRIX
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'tilt=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # tilt angle
+    if test:
+        tilt = float(test.group(1))
+        if not (tilt > -180) or not (tilt < 180):
+            errMsg(
+                "TILT=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
+            ) % tilt
+            sys.exit(1)
+        tiltang = numpy.multiply(numpy.py, numpy.divide(tilt, 180))
+        sintilt = numpy.sin(tiltang)
+        costilt = numpy.cos(tiltang)
+        Rt0 = [1, 0, 0]
+        Rt1 = [0, costilt, sintilt]
+        Rt2 = [0, -sintilt, costilt]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rt0, Rt1, Rt2, R0, R1, R2)
+        #again resets the matrix ... hmm
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'roll=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # roll angle
+    if test:
+        roll = float(test.group(1))
+        if not (roll > -180) or not (roll < 180):
+            errMsg(
+                "ROLL=%s MUST BE LESS THAN 180 AND GREATER THAN -180"
+            ) % roll
+
+        rollang = numpy.multiply(numpy.py, numpy.divide(roll, 180))
+        sinroll = numpy.sin(rollang)
+        cosroll = numpy.cos(rollang)
+        Rr0 = [cosroll, sinroll, 0]
+        Rr1 = [-sinroll, cosroll, 0]
+        Rr2 = [0, 0, 1]
+        # do matrix multiply to get new rotation matrix
+        newmat = matmul3(Rr0, Rr1, Rr2, R0, R1, R2)
+        R0 = newmat[0]
+        R1 = newmat[1]
+        R2 = newmat[2]
+    test = re.match(r'pef=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # pef angle
+    if test:
+        pef = float(test.group(1))
+        if not (pef > 0) or not (pef < numpy.pi):
+            errMsg(
+                "PEF=%s MUST BE LESS THAN PI AND GREATER THAN 0"
+            ) % pef
+    test = re.match(r'idx=(\d+\.\d+)', arg, flags=re.IGNORECASE)  # idx angle
+    if test:
+        idx = float(test.group(1))
+        if not (idx >= 0):
+            errMsg(
+                "IDX=%s MUST BE 0 or GREATER"
+            ) % idx
+        di = idx
+    test = re.match(r'idy=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        idy = float(test.group(1))
+        if not (idy >= 0):
+            errMsg(
+                "IDY=%s MUST BE 0 or GREATER"
+            ) % idy
+        dj = idy
+    test = re.match(r'odx=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        odx = float(test.group(1))
+        if not (odx >= 0):
+            errMsg(
+                "ODX=%s MUST BE 0 or GREATER"
+            ) % odx
+        du = odx
+    test = re.match(r'ody=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        ody = float(test.group(1))
+        if not (ody > 0):
+            errMsg(
+                "ODY=%s MUST BE 0 or GREATER"
+            ) % ody
+        dv = ody
+    test = re.match(r'zoom=(\d+\.\d+)', arg, flags=re.IGNORECASE)
+    if test:
+        zoom = float(test.group(1))
+        if not (zoom < -1 or zoom > 1):
+            errMsg(
+                "ODY=%s MUST BE GREATER THAN 1 OR LESS THAN -1"
+            ) % zoom
+    test = re.match(r'bgcolor=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        bgcolor = test.group(1)
+    test = re.match(r'skycolor=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        skycolor = test.group(1)
+    test = re.match(r'vp=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        vp = test.group(1).lower()
+        if vp not in [
+                "background", "dither", "edge",
+                "mirror", "random", "tile", "transparent"
+        ]:
+            errMsg("VP=%s IS NOT A VALID VALUE" % (vp))
+    test = re.match(r'auto=([^ ]+)', arg, flags=re.IGNORECASE)
+    if test:
+        auto = test.group(1).lower()
+        if auto not in [
+                "c", "zc", "out"
+        ]:
+            errMsg("AUTO=%s IS NOT A VALID VALUE" % (auto))
+#
+# get infile and outfile
+infile = sys.argv[-2]
+outfile = sys.argv[-1]
+for a in [infile, outfile]:
+    if re.match(r'.*=.*', a):
+        errMsg("%s is not a valid filepath" % a)
+
+
+# setup temporary images and auto delete upon exit
+# use mpc/cache to hold input image temporarily in memory
+pid = os.getpid()
+
+tmpA = "%s/3Drotate_%s.mpc" % (directory, pid)
+tmpB = "%s/3Drotate_%s.cache" % (directory, pid)
+os.unlink(tmpA, tmpB)
+os.unlink(tmpA, tmpB)
+
+# test that infile provided
+if not os.path.exists(infile):
+    errMsg("NO INPUT FILE SPECIFIED")
+# test that outfile provided
+if not outfile:
+    errMsg("NO OUTPUT FILE SPECIFIED")
+
+cmd = ['convert', '-quiet', '-regard-warnings', infile, '+repage', tmpA]
+call_cmd(cmd)
+if not pef:
+    pef = 1
+
+# get input image width and height
+width, height = imagesize()
+maxwidth = width - 1
+maxheight = height - 1
+
+# deal with auto adjustments to values
+if auto == "zc":
+    du = 0
+    dv = 0
+    zoom = 1
+elif auto == "c":
+    du = 0
+    dv = 0
+
+# convert offsets of rotation point to relative to pixel 0,0
+di = numpy.divide(numpy.add(di, numpy.divide(numpy.subtract(width, 1), 2)), 1)
+dj = numpy.divide(
+        numpy.add(dj, numpy.divide(numpy.subtract(height, 1), 2)), 1)
+du = numpy.divide(du, 1)
+dv = numpy.divide(dv, 1)
+
+# convert zoom to scale factors
+if zoom >= 1:
+    sx = numpy.divide(1, zoom)
+    sy = sx
+elif zoom <= -1:
+    sx = numpy.divide(-zoom, 1)
+    sy = sx
+
+#{{{explanation
+# Consider the picture placed on the Z=0 plane and the camera a distance
+# Zc=f above the picture plane looking straight down at the image center.
+# Now the perspective equations (in 3-D) are defined as (x,y,f) = M (X',Y',Z'),
+# where the camera orientation matrix M is the identity matrix but with M22=-1
+# because the camera is looking straight down along -Z.
+# Thus a reflection transformation relative to the ground plane coordinates.
+# Let the camera position Zc=f=(sqrt(ins*ins + inl*inl)) / ( 2 tan(fov/2) )
+# Now we want to rotate the ground points corresponding to the picture corners.
+# The basic rotation is (X',Y',Z') = R (X,Y,0), where R is the rotation matrix
+# involving pan, tilt and roll.
+# But we need to convert (X,Y,0) to (X,Y,1) and also to offset for Zc=f
+# First we note that (X,Y,0) = (X,Y,1) - (0,0,1)
+# Thus the equation becomes
+# (x,y,f) = M {R [(X,Y,1) - (0,0,1)] - (0,0,Zc)} = MT (X,Y,1)
+# But R [(X,Y,1) - (0,0,1)] = R [II (X,Y,1) - S (X,Y,1)] = R (II-S) (X,Y,1),
+# where II is the identity matrix and S is an all zero matrix except for S22=1.
+# Thus (II-S) is the identity matrix with I22=0 and
+# RR = R (II-S) is just R with the third column all zeros.
+# Thus we get (x,y,f) = M {RR (X,Y,1) - (0,0,Zc)}.
+# But M {RR (X,Y,1) - (0,0,Zc)} = M {RR(X,Y,1) - D (X,Y,1)}, where
+# D is an all zero matrix with D22 = Zc = f.
+# So that we get M (RR-D) (X,Y,1) = MT (X,Y,1), where
+# where T is just R with the third column (0,0,-f), i.e. T02=0, T12=0, T22=-f
+# But we need to allow for scaling and offset of the output coordinates and
+# conversion from (x,y,f) to (u,v,1)=O and conversion of input coordinates
+# from (X,Y,1) to (i,j,1)=I.
+# Thus the forward transformation becomes AO=MTBI or O=A'MTBI or O=PI,
+# where prime means inverse.
+# However, to do the scaling of the output correctly, need to offset by the
+# input plus output offsets, then scale, which is all put into A'.
+# Thus the forward transformation becomes AO=MTBI or O=A'MTBI where A'=Ai
+# but we will merge A'M into Aim
+# Thus the inverse transform becomes
+# I=QO where Q=P'
+# A=output scaling, offset and conversion matrix
+# B=input offset and conversion matrix
+# (scaling only needs to be done in one place)
+# M=camera orientation matrix
+# R=image rotation matrix Rroll Rtilt Rpan
+# T=matrix that is R but R33 offset by f + 1
+# O=output coords vector (i,j,1)
+# I=input coords vector (u,v,1)=(is,il,1)
+# P=forward perspective transformation matrix
+# Q=inverse perspective transformation matrix
+#
+# For a 35 mm camera whose film format is 36mm wide
+# and 24mm tall, when the focal length
+# is equal to the diagonal, the field of view is 53.13 degrees and this is
+# considered a normal view equivalent to the human eye.
+# See http://www.panoramafactory.com/equiv35/equiv35.html
+# Max limit on dfov is 180 degrees (pef=3.19)
+# where get single line like looking at picture on edge.
+# Above this limit the picture becomes like the angles get reversed.
+# Min limit on dfov seems to be slightly greater than zero degrees.
+# Practical limits on dfov depend upon orientation angles.
+# For tilt=45, this is about 2.5 dfov (pef=2.5).
+# Above this, some parts of the picture
+# that are cut off at the bottom, get wrapped and stretched in the 'sky'.
+#}}}
+
+dfov = numpy.divide(
+        numpy.multpily(180, numpy.arctan(36/24)), numpy.pi)
+pfact = 1
+if pef == 0:  # FIXME
+    pfact = numpy.divide(0.01, dfov)
+else:
+    pfact = pef
+
+#compute new field of view based upon pef (pfact)
+dfov = numpy.multiply(pfact, dfov)
+dfov2 = numpy.divide(dfov, 2)
+arg = numpy.multiply(numpy.pi, numpy.divide(dfov, 180))
+sfov = numpy.sin(arg)
+cfov = numpy.cos(arg)
+tfov = numpy.divide(sfov, cfov)
+
+
+# calculate focal length in same units as wall (picture) using dfov
+diag = numpy.sqrt(
+        numpy.multiply(width, width), numpy.multiply(height, height))
+focal = numpy.divide(diag, numpy.multiply(2, tfov))
+
+# calculate forward transform matrix Q
+
+# define the input offset and conversion matrix
+dim = -di
+B0 = [1, 0, dim]
+B1 = [0, -1, dj]
+B2 = [0, 0, 1]
+
+# define the output scaling, offset and conversion
+# matrix inverse Ai and merge with M
+# to become Aim
+#A0=($sx 0 $sx*(-$du-$di))
+#A1=(0 -$sy $sy*($dv+$dj))
+#A2=(0 0 -$focal)
+#M0=(1 0 0)
+#M1=(0 1 0)
+#M2=(0 0 -1)
+aim00 = numpy.divide(1, sx)
+aim02 = -numpy.divide(
+        numpy.multiply(sx, numpy.add(di, du)), numpy.multiply(sx, focal))
+aim11 = numpy.divide(-1, sy)
+aim12 = -numpy.divide(
+        numpy.multiply(sy, numpy.add(dj, dv)), numpy.multiply(sy, focal))
+aim22 = -numpy.divide(1, focal)
+Aim0 = [aim00, 0, aim02]
+Aim1 = [0, aim11, aim12]
+Aim2 = [0, 0, aim22]
+
+# now do successive matrix multiplies
+# from right towards left of main equation P=A'RB
+
+# convert R to T by setting T02=T12=0 and T22=-f
+T0 = [R0[0], R0[1], 0]
+T1 = [R1[0], R1[1], 0]
+T2 = [R2[0], R2[1], -focal]
+
+# multiply T x B = P
+p_mat = matmul3(T0, T1, T2, B0, B1, B2)
+
+# multiply Aim x P = P
+newmat = matmul3(Aim0, Aim1, Aim2, p_mat[0], p_mat[1], p_mat[2])
+
+#    # the resulting P matrix is now
+#    the perspective coefficients for the inverse transformation
+#    P00= newmat[0][0]
+#    P01= newmat[0][1]
+#    P02= newmat[0][2]
+#    P10= newmat[1][0]
+#    P11= newmat[1][1]
+#    P12= newmat[1][2]
+#    P20= newmat[2][0]
+#    P21= newmat[2][1]
+#    P22= newmat[2][2]
+
+# project input corners to output domain
+#echo "UL"
+i = 0
+j = 0
+u1, v1 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
+
+i = maxwidth
+j = 0
+u2, v2 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
+
+i = maxwidth
+j = maxheight
+
+u3, v3 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
+i = 0
+j = maxheight
+
+u4, v4 = forwardProject(newmat, i, j)  # using Aim matrix and p_mat
+
+# deal with adjustments for auto settings
+# first get the bounding box dimensions
+uArr = [u1, u2, u3, u4]
+vArr = [v1, v2, v3, v4]
+umin = 1000000
+umax = -1000000
+vmin = 1000000
+vmax = -1000000
+for index in range(0, 4):
+    if uArr[index] < umin:
+        umin = uArr[index]
+    if uArr[index] > umax:
+        umax = uArr[index]
+    if vArr[index] < vmin:
+        vmin = vArr[index]
+    if vArr[index] > vmax:
+        vmax = vArr[index]
+delu = numpy.add(numpy.subtract(umax, umin), 1)
+delv = numpy.add(numpy.subtract(vmax, vmin), 1)
+if auto == "c":
+    offsetu = numpy.divide(numpy.subtract(width, delu), 2)
+    offsetv = numpy.divide(numpy.subtract(height, delv), 2)
+    #FIXME may need to cast as int variables below (u1-v4)
+    u1 = numpy.add(offsetu, numpy.subtract(u1, umin))
+    v1 = numpy.add(offsetv, numpy.subtract(v1, vmin))
+    u2 = numpy.add(offsetu, numpy.subtract(u2, umin))
+    v2 = numpy.add(offsetv, numpy.subtract(v2, vmin))
+    u3 = numpy.add(offsetu, numpy.subtract(u3, umin))
+    v3 = numpy.add(offsetv, numpy.subtract(v3, vmin))
+    u4 = numpy.add(offsetu, numpy.subtract(u4, umin))
+    v4 = numpy.add(offsetv, numpy.subtract(v4, vmin))
+elif auto == "zc":
+    if delu > delv:
+        _del = delu
+        offsetu = 0
+        offsetv = numpy.divide(
+            numpy.subtract(
+                height, numpy.divide(numpy.multpily(delv, width), delu)), 2)
+    else:
+        _del = delv
+        offsetu = numpy.divide(
+            numpy.subtract(
+                width, numpy.divide(numpy.multpily(delu, height), delv)), 2)
+        offsetv = 0
+u1 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u1, umin), numpy.divide(width, _del)))
+v1 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v1, vmin), numpy.divide(height, _del)))
+u2 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u2, umin), numpy.divide(width, _del)))
+v2 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v2, vmin), numpy.divide(height, _del)))
+u3 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u3, umin), numpy.divide(width, _del)))
+v3 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v3, vmin), numpy.divide(height, _del)))
+u4 = numpy.add(
+    offsetu,
+    numpy.multpily(numpy.subtract(u4, umin), numpy.divide(width, _del)))
+v4 = numpy.add(
+    offsetv,
+    numpy.multpily(numpy.subtract(v4, vmin), numpy.divide(height, _del)))
+#
+# now do the perspective distort
+if auto == "out":
+    distort = "+distort"
+else:
+    distort = "-distort"
+
+#im_version=`convert -list configure | \
+#	sed '/^LIB_VERSION_NUMBER /!d; s//,/;  s/,/,0/g;  s/,0*\([0-9][0-9]\)/\1/g'\
+#	| head -n 1`
+#if [ "$im_version" -lt "06030600" ]
+#	then
+#	convert $tmpA -virtual-pixel $vp -background $bgcolor \
+#	-mattecolor $skycolor $distort Perspective \
+#	"0,0 $maxwidth,0 $maxwidth,$maxheight 0,$maxheight  \
+#	$u1,$v1 $u2,$v2 $u3,$v3 $u4,$v4" $outfile
+#else
+cmd = [
+    "convert", tmpA, "-virtual-pixel", vp, "-background", bgcolor,
+    "-mattecolor", skycolor, distort, "Perspective",
+    "0,0", "%s,%s" % (u1, v1),
+    "%s,0" % (maxwidth), "%s,%s" % (u2, v2),
+    "%s,%s" % (maxwidth, maxheight),
+    "%s,%s" % (u3, v3),
+    "%s,%s" % (0, maxheight),
+    "%s,%s" % (u4, v4),
+    outfile
+]
+call_cmd(cmd)