almost done

1 files changed, 230 insertions, 265 deletions

diff --git a/bin/3Drotate.py b/bin/3Drotate.py
index 63281b7..962a285 100755
--- a/bin/3Drotate.py
+++ b/bin/3Drotate.py
@@ -261,7 +261,7 @@ else:
             sys.exit(1)
         if arg == "-":
             break  #fixme
-        test = re.match(r'pan=(\d+\.\d+)', arg)
+        test = re.match(r'pan=(\d+\.\d+)', arg, flags=re.IGNORECASE)
         if test:
             pan = float(test.group(1))
             if (pan > 180) or (pan < -180):
@@ -277,10 +277,10 @@ else:
             # 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) # tilt angle
+            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):
@@ -299,282 +299,247 @@ else:
             R0 = newmat[0]
             R1 = newmat[1]
             R2 = newmat[2]
-		  roll[=]*)    # roll angle
-					   arg="$1="
-					   roll=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   roll=`echo "$roll" | sed 's/^[+]\(.*\)$/\1/'`
-					   # rolltest>0 if floating point number; otherwise rolltest=0
-					   testFloat "$roll"; rolltest=$floatresult
-					   rolltestA=`echo "$roll < - 180" | bc`
-					   rolltestB=`echo "$roll > 180" | bc`
-					   [ $rolltest -eq 0 ] && errMsg "roll=$roll IS NOT A NUMBER"
-					   [ $rolltestA -eq 1 -o $rolltestB -eq 1 ] && errMsg "ROLL=$roll MUST BE GREATER THAN -180 AND LESS THAN +180"
-					   rollang=`echo "scale=10; $pi * $roll / 180" | bc`
-					   sinroll=`echo "scale=10; s($rollang)" | bc -l`
-					   sinrollm=`echo "scale=10; - $sinroll" | bc`
-					   cosroll=`echo "scale=10; c($rollang)" | bc -l`
-					   Rr0=($cosroll $sinroll 0)
-					   Rr1=($sinrollm $cosroll 0)
-					   Rr2=(0 0 1)
-					   # do matrix multiply to get new rotation matrix
-					   matmul3 "${Rr0[*]}" "${Rr1[*]}" "${Rr2[*]}" "${R0[*]}" "${R1[*]}" "${R2[*]}"
-					   R0=(${P0[*]})
-					   R1=(${P1[*]})
-					   R2=(${P2[*]})
-					   ;;
-		   pef[=]*)    # pef
-					   arg="$1="
-					   pef=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   pef=`echo "$pef" | sed 's/^[+]\(.*\)$/\1/'`
-					   # peftest>0 if floating point number; otherwise peftest=0
-					   testFloat "$pef"; peftest=$floatresult
-					   peftestA=`echo "$pef < 0" | bc`
-					   peftestB=`echo "$pef > 3.19" | bc`
-					   [ $peftest -eq 0 ] && errMsg "PEF=$pef IS NOT A NUMBER"
-					   ;;
-	       idx[=]*)    # input x shift
-					   arg="$1="
-					   di=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   di=`echo "$di" | sed 's/^[+]\(.*\)$/\1/'`
-					   # ditest>0 if floating point number; otherwise ditest=0
-					   testFloat "$di"; ditest=$floatresult
-					   [ $ditest -eq 0 ] && errMsg "ISHIFTX=$di IS NOT A NUMBER"
-					   ;;
-	       idy[=]*)    # input y shift
-					   arg="$1="
-					   dj=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   dj=`echo "$dj" | sed 's/^[+]\(.*\)$/\1/'`
-					   # djtest>0 if floating point number; otherwise ditest=0
-					   testFloat "$dj"; djtest=$floatresult
-					   [ $djtest -eq 0 ] && errMsg "ISHIFTY=$dj IS NOT A NUMBER"
-					   ;;
-	       odx[=]*)    # output x shift
-					   arg="$1="
-					   du=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   du=`echo "$du" | sed 's/^[+]\(.*\)$/\1/'`
-					   # dutest>0 if floating point number; otherwise ditest=0
-					   testFloat "$du"; dutest=$floatresult
-					   [ $dutest -eq 0 ] && errMsg "OSHIFTX=$du IS NOT A NUMBER"
-					   ;;
-	       ody[=]*)    # output y shift
-					   arg="$1="
-					   dv=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   dv=`echo "$dv" | sed 's/^[+]\(.*\)$/\1/'`
-					   # dvtest>0 if floating point number; otherwise ditest=0
-					   testFloat "$dv"; dvtest=$floatresult
-					   [ $dvtest -eq 0 ] && errMsg "OSHIFTY=$dv IS NOT A NUMBER"
-					   ;;
-		  zoom[=]*)    # output zoom
-					   arg="$1="
-					   zoom=`echo "$arg" | cut -d= -f2`
-					   # function bc does not seem to like numbers starting with + sign, so strip off
-					   zoom=`echo "$zoom" | sed 's/^[+]\(.*\)$/\1/'`
-					   # zoomtest>0 if floating point number; otherwise peftest=0
-					   testFloat "$zoom"; zoomtest=$floatresult
-					   zoomtest=`echo "$zoom < 1 && $zoom > -1" | bc`
-					   [ $zoomtest -eq 1 ] && errMsg "ZOOM=$zoom MUST BE GREATER THAN 1 OR LESS THAN -1"
-					   ;;
-	   bgcolor[=]*)    # output background color
-					   arg="$1="
-					   bgcolor=`echo "$arg" | cut -d= -f2`
-					   ;;
-	  skycolor[=]*)    # output sky color
-					   arg="$1="
-					   skycolor=`echo "$arg" | cut -d= -f2`
-					   ;;
-	        vp[=]*)    # virtual pixel method
-					   arg="$1="
-					   vp=`echo "$arg" | cut -d= -f2`
-					   [ "$vp" != "background" -a "$vp" != "dither" -a "$vp" != "edge" -a "$vp" != "mirror" -a "$vp" != "random" -a "$vp" != "tile" -a "$vp" != "transparent" ] && errMsg "VP=$vp IS NOT A VALID VALUE"
-					   ;;
-	   	  auto[=]*)    # output background color
-					   arg="$1="
-					   auto=`echo "$arg" | cut -d= -f2`
-					   [ "$auto" != "c" -a "$auto" != "zc" -a "$auto" != "out" ] && errMsg "AUTO=$auto IS NOT A VALID VALUE"
-					   ;;
-		     *[=]*)    # not valid
-					   errMsg "$1 IS NOT A VALID ARGUMENT"
-					   ;;
-		     	 *)    # end of arguments
-					   break
-					   ;;
-			esac
-			shift   # next option
-	done
-	#
+        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=$1
-	outfile=$2
-fi
+    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)
+    # setup temporary images and auto delete upon exit
+    # use mpc/cache to hold input image temporarily in memory
+    pid = os.getpid()
 
-# test that infile provided
-[ "$infile" = "" ] && errMsg "NO INPUT FILE SPECIFIED"
-# test that outfile provided
-[ "$outfile" = "" ] && errMsg "NO OUTPUT FILE SPECIFIED"
+    tmpA = "%s/3Drotate_%s.mpc" % (directory, pid)
+    tmpB = "%s/3Drotate_%s.cache" % (directory, pid)
+    os.unlink(tmpA,tmpB)
+    os.unlink(tmpA,tmpB)
 
-if convert -quiet -regard-warnings "$infile" +repage "$tmpA"
-	then
-	[ "$pef" = "" ] && pef=1
-else
-	errMsg "--- FILE $infile DOES NOT EXIST OR IS NOT AN ORDINARY FILE, NOT READABLE OR HAS ZERO SIZE ---"
-fi
+    # 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")
 
-# get input image width and height
-imagesize
-maxwidth=`expr $width - 1`
-maxheight=`expr $height - 1`
+    cmd = ['convert', '-quiet', '-regard-warnings', infile, '+repage', tmpA]
+    call_cmd(cmd)
+    if not pef:
+        pef = 1
 
-# deal with auto adjustments to values
-if [ "$auto" = "zc" ]
-	then
-	du=0
-	dv=0
-	zoom=1
-elif [ "$...auto" = "c" ]
-	then
-	du=0
-	dv=0
-fi
+    # get input image width and height
+    imagesize()
+    maxwidth = width - 1
+    maxheight = height - 1
 
-# convert offsets of rotation point to relative to pixel 0,0
-di=`echo "scale=10; ($di + (($width - 1) / 2)) / 1" | bc`
-dj=`echo "scale=10; ($dj + (($height - 1) / 2)) / 1" | bc`
-du=`echo "scale=10; $du / 1" | bc`
-dv=`echo "scale=10; $dv / 1" | bc`
+    # deal with auto adjustments to values
+    if auto == "zc":
+        du = 0
+        dv = 0
+        zoom = 1
+    elif auto == "c":
+        du = 0
+        dv = 0
 
-# convert zoom to scale factors
-if [ `echo "$zoom >= 1" | bc` -eq 1 ]
-	then
-	sx=`echo "scale=10; 1 / $zoom" | bc`
-	sy=$sx
-elif [ `echo "$zoom <= -1" | bc` -eq 1 ]
-	then
-	sx=`echo "scale=10; - $zoom / 1" | bc`
-	sy=$sx
-fi
+    # 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)
 
-#{{{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'.
-#}}}
+    # 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
 
-dfov=`echo "scale=10; 180 * a(36/24) / $pi" | bc -l`
-if [ "$pef" = "" ]
-	then
-	pfact=1
-elif [ "$pef" = "0" ]
-	then
-	pfact=`echo "scale=10; 0.01 / $dfov" | bc`
-else
-	pfact=$pef
-fi
-#maxpef=`echo "scale=5; 180 / $dfov" | bc`
-#echo "maxpef=$maxpef"
+    #{{{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)
 
-#compute new field of view based upon pef (pfact)
-dfov=`echo "scale=10; $pfact * $dfov" | bc`
-dfov2=`echo "scale=10; $dfov / 2" | bc`
-arg=`echo "scale=10; $pi * $dfov2 / 180" | bc`
-sfov=`echo "scale=10; s($arg)" | bc -l`
-cfov=`echo "scale=10; c($arg)" | bc -l`
-tfov=`echo "scale=10; $sfov / $cfov" | bc -l`
-#echo "tfov=$tfov"
 
-# calculate focal length in same units as wall (picture) using dfov
-diag=`echo "scale=10; sqrt(($width * $width) + ($height * $height))" | bc`
-focal=`echo "scale=10; ($diag / (2 * $tfov))" | bc -l`
-#echo "focal=$focal"
+    # 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=`echo "scale=10; - $di" | bc`
-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=`echo "scale=10; 1 / $sx" | bc`
-aim02=`echo "scale=10; -($sx * ($di + $du)) / ($sx * $focal)" | bc`
-aim11=`echo "scale=10; -1 / $sy" | bc`
-aim12=`echo "scale=10; -($sy * ($dj + $dv)) / ($sy * $focal)" | bc`
-aim22=`echo "scale=10; -1 / $focal" | bc`
-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= -1 / $sy" | bc`
+    aim12= -($sy * ($dj + $dv)) / ($sy * $focal)" | bc`
+    aim22= -1 / $focal" | bc`
+    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