From 95cf222e29dfce0f2c3e2771912aa007c64cd874 Mon Sep 17 00:00:00 2001 From: yo mama Date: Sat, 24 Jun 2017 20:21:21 -0700 Subject: added 3Drotate script --- bin/3Drotate | 39 ---- bin/3Drotate.py | 696 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ bin/grid | 9 +- 3 files changed, 701 insertions(+), 43 deletions(-) create mode 100755 bin/3Drotate.py (limited to 'bin') diff --git a/bin/3Drotate b/bin/3Drotate index 227b92d..55b50e6 100755 --- a/bin/3Drotate +++ b/bin/3Drotate @@ -757,45 +757,6 @@ v4=$vv #u5=$uu #v5=$vv -# unused -: ' -# Now invert P to get Q for the inverse perspective transformation -# Use the Method of the Adjoint Matrix = transpose of matrix of cofactors divided by the determinant -# M3inverse $P00 $P01 $P02 $P10 $P11 $P12 $P20 $P21 $P22 -# -# project output corners to input domain -# UL -#echo "UL 0,0" -#u=$u1 -#v=$v1 -#echo "u,v=$u,$v" -#inverseProject $u $v -#echo "i,j=$ii,$jj" -#echo "UR 255,0" -#u=$u2 -#v=$v2 -#echo "u,v=$u,$v" -#inverseProject $u $v -#echo "i,j=$ii,$jj" -#echo "BR 255,255" -#u=$u3 -#v=$v3 -#echo "u,v=$u,$v" -#inverseProject $u $v -#echo "i,j=$ii,$jj" -#echo "BL 0,255" -#u=$u4 -#v=$v4 -#echo "u,v=$u,$v" -#inverseProject $u $v -#echo "i,j=$ii,$jj" -#echo "C 127.5,127.5" -#u=$u5 -#v=$v5 -#echo "u,v=$u,$v" -#inverseProject $u $v -#echo "i,j=$ii,$jj" -' # deal with adjustments for auto settings # first get the bounding box dimensions diff --git a/bin/3Drotate.py b/bin/3Drotate.py new file mode 100755 index 0000000..63281b7 --- /dev/null +++ b/bin/3Drotate.py @@ -0,0 +1,696 @@ +#!/usr/bin/python2.7 +import getopt +import numpy +import sys +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 +dir = "." # 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 forwardProjrect(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 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(" ") + ps = 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) +elif (len_args > 15): + errMsg("--- TOO MANY ARGUMENTS WERE PROVIDED ---") +else: + for arg in sys.argv: + if arg == "-h": + usage2() + sys.exit(1) + if arg == "-": + break #fixme + test = re.match(r'pan=(\d+\.\d+)', arg) + 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) # 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] + 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 + # + # get infile and outfile + infile=$1 + outfile=$2 +fi + +# 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 +[ "$infile" = "" ] && errMsg "NO INPUT FILE SPECIFIED" +# test that outfile provided +[ "$outfile" = "" ] && errMsg "NO OUTPUT FILE SPECIFIED" + +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 + +# get input image width and height +imagesize +maxwidth=`expr $width - 1` +maxheight=`expr $height - 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 + +# 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` + +# 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 + +#{{{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=`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" + +#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 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 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) + +# 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 +focalm=`echo "scale=10; - $focal" | bc` +T0=(${R0[0]} ${R0[1]} 0) +T1=(${R1[0]} ${R1[1]} 0) +T2=(${R2[0]} ${R2[1]} $focalm) + +# 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) +index=0 +umin=1000000 +umax=-1000000 +vmin=1000000 +vmax=-1000000 +while [ $index -lt 4 ] + do + [ `echo "${uArr[$index]} < $umin" | bc` -eq 1 ] && umin=${uArr[$index]} + [ `echo "${uArr[$index]} > $umax" | bc` -eq 1 ] && umax=${uArr[$index]} + [ `echo "${vArr[$index]} < $vmin" | bc` -eq 1 ] && vmin=${vArr[$index]} + [ `echo "${vArr[$index]} > $vmax" | bc` -eq 1 ] && vmax=${vArr[$index]} + index=`expr $index + 1` +done +delu=`echo "scale=10; $umax - $umin + 1" | bc` +delv=`echo "scale=10; $vmax - $vmin + 1" | bc` +if [ "$auto" = "c" ] + then + 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=`echo "scale=0; $offsetu + (($u1 - $umin) * $width / $del)" | bc` + v1=`echo "scale=0; $offsetv + (($v1 - $vmin) * $height / $del)" | bc` + u2=`echo "scale=0; $offsetu + (($u2 - $umin) * $width / $del)" | bc` + v2=`echo "scale=0; $offsetv + (($v2 - $vmin) * $height / $del)" | bc` + u3=`echo "scale=0; $offsetu + (($u3 - $umin) * $width / $del)" | bc` + v3=`echo "scale=0; $offsetv + (($v3 - $vmin) * $height / $del)" | bc` + u4=`echo "scale=0; $offsetu + (($u4 - $umin) * $width / $del)" | bc` + v4=`echo "scale=0; $offsetv + (($v4 - $vmin) * $height / $del)" | bc` +fi +# +# now do the perspective distort +if [ "$auto" = "out" ] + then + distort="+distort" +else + distort="-distort" +fi + +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 diff --git a/bin/grid b/bin/grid index 01678ac..43762a2 100755 --- a/bin/grid +++ b/bin/grid @@ -186,13 +186,14 @@ fi # use mpc/cache to hold input image temporarily in memory tmpA="$dir/profile_$$.mpc" tmpB="$dir/profile_$$.cache" -trap "rm -f $tmpA $tmpB; exit 0" 0 -trap "rm -f $tmpA $tmpB; exit 1" 1 2 3 15 +#trap "rm -f $tmpA $tmpB; exit 0" 0 +#trap "rm -f $tmpA $tmpB; exit 1" 1 2 3 15 # -if convert -quiet -regard-warnings "$infile" +repage "$tmpA" - then +convert "$infile" +repage "$tmpA" +if [ $? -eq 0 ]; + then : 'do nothing - continue processing below' else errMsg "--- FILE $infile DOES NOT EXIST OR IS NOT AN ORDINARY FILE, NOT READABLE OR HAS ZERO SIZE ---" -- cgit v1.2.3-70-g09d2 From 91ca7d68e4caa9d6875549bfd2c329b9d0ef48ec Mon Sep 17 00:00:00 2001 From: yo mama Date: Sun, 25 Jun 2017 16:35:58 -0700 Subject: almost done --- bin/3Drotate.py | 525 ++++++++++++++++++++++++++------------------------------ 1 file changed, 245 insertions(+), 280 deletions(-) (limited to 'bin') 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 - -# 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 -[ "$infile" = "" ] && errMsg "NO INPUT FILE SPECIFIED" -# test that outfile provided -[ "$outfile" = "" ] && errMsg "NO OUTPUT FILE SPECIFIED" - -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 - -# get input image width and height -imagesize -maxwidth=`expr $width - 1` -maxheight=`expr $height - 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 - -# 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` - -# 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 - -#{{{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=`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" - -#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 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 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) + 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 + 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= -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 -- cgit v1.2.3-70-g09d2 From ea3138b1f6cad5d1958ca86030846939f926a267 Mon Sep 17 00:00:00 2001 From: yo mama Date: Sun, 25 Jun 2017 17:11:57 -0700 Subject: so close --- bin/3Drotate.py | 197 ++++++++++++++++++++++++++++++-------------------------- 1 file changed, 107 insertions(+), 90 deletions(-) (limited to 'bin') diff --git a/bin/3Drotate.py b/bin/3Drotate.py index 962a285..7669ffa 100755 --- a/bin/3Drotate.py +++ b/bin/3Drotate.py @@ -534,88 +534,89 @@ else: 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 - -# convert R to T by setting T02=T12=0 and T22=-f -focalm=`echo "scale=10; - $focal" | bc` -T0=(${R0[0]} ${R0[1]} 0) -T1=(${R1[0]} ${R1[1]} 0) -T2=(${R2[0]} ${R2[1]} $focalm) - -# 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) -index=0 -umin=1000000 -umax=-1000000 -vmin=1000000 -vmax=-1000000 -while [ $index -lt 4 ] - do - [ `echo "${uArr[$index]} < $umin" | bc` -eq 1 ] && umin=${uArr[$index]} - [ `echo "${uArr[$index]} > $umax" | bc` -eq 1 ] && umax=${uArr[$index]} - [ `echo "${vArr[$index]} < $vmin" | bc` -eq 1 ] && vmin=${vArr[$index]} - [ `echo "${vArr[$index]} > $vmax" | bc` -eq 1 ] && vmax=${vArr[$index]} - index=`expr $index + 1` -done -delu=`echo "scale=10; $umax - $umin + 1" | bc` -delv=`echo "scale=10; $vmax - $vmin + 1" | bc` -if [ "$auto" = "c" ] - then - 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` + 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 = `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 ] @@ -628,14 +629,30 @@ elif [ "$auto" = "zc" ] offsetu=`echo "scale=10; ($width - ($delu * $height / $delv)) / 2" | bc` offsetv=0 fi - u1=`echo "scale=0; $offsetu + (($u1 - $umin) * $width / $del)" | bc` - v1=`echo "scale=0; $offsetv + (($v1 - $vmin) * $height / $del)" | bc` - u2=`echo "scale=0; $offsetu + (($u2 - $umin) * $width / $del)" | bc` - v2=`echo "scale=0; $offsetv + (($v2 - $vmin) * $height / $del)" | bc` - u3=`echo "scale=0; $offsetu + (($u3 - $umin) * $width / $del)" | bc` - v3=`echo "scale=0; $offsetv + (($v3 - $vmin) * $height / $del)" | bc` - u4=`echo "scale=0; $offsetu + (($u4 - $umin) * $width / $del)" | bc` - v4=`echo "scale=0; $offsetv + (($v4 - $vmin) * $height / $del)" | bc` + 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 # # now do the perspective distort -- cgit v1.2.3-70-g09d2 From a69680d9347372ffadbb52e455311f6d20ccaf57 Mon Sep 17 00:00:00 2001 From: yo mama Date: Sun, 25 Jun 2017 17:48:11 -0700 Subject: new --- bin/3Drotate.py | 896 +++++++++++++++++++++++++++------------------------- bin/3Drotate.py.bak | 716 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 1183 insertions(+), 429 deletions(-) create mode 100755 bin/3Drotate.py.bak (limited to 'bin') 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 - - 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 - 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) + 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: - 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 = `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 + _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 - -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 +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) 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) -- cgit v1.2.3-70-g09d2 From abff780ddb60a5218c3a463a7e257024e02476a2 Mon Sep 17 00:00:00 2001 From: yo mama Date: Sun, 25 Jun 2017 18:06:08 -0700 Subject: awesome, debugging --- bin/3Drotate.py | 81 ++++++++++++++++++++++++++++++++------------------------- 1 file changed, 45 insertions(+), 36 deletions(-) (limited to 'bin') diff --git a/bin/3Drotate.py b/bin/3Drotate.py index a94a7a9..28f1443 100755 --- a/bin/3Drotate.py +++ b/bin/3Drotate.py @@ -156,6 +156,7 @@ foolproof. Use At Your Own Risk. #{{{ defaults # set default value # rotation angles and rotation matrix +auto = None pan = 0 tilt = 0 roll = 0 @@ -220,7 +221,7 @@ def forwardProject(mat, ii, jj): def errMsg(s): - sys.stderr.write("%s\n") + sys.stderr.write("%s\n" % s) sys.exit(1) @@ -246,8 +247,7 @@ def M3inverse(mat): return numpy.linalg.inv(mat) -def call_cmd(s): - cmd = s.split(" ") +def call_cmd(cmd): process = Popen(cmd, stdout=PIPE, stderr=PIPE) out, err = process.communicate() errcode = process.returncode @@ -258,8 +258,8 @@ def call_cmd(s): # 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))) + width = int(call_cmd(["identify", "-format", "%w", "%s" % (filepath)])) + height = int(call_cmd(["identify", "-format", "%h", "%s" % (filepath)])) return width, height @@ -417,8 +417,9 @@ 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) +for a in [tmpA, tmpB]: + if os.path.exists(a): + os.unlink(a) # test that infile provided if not os.path.exists(infile): @@ -433,7 +434,7 @@ if not pef: pef = 1 # get input image width and height -width, height = imagesize() +width, height = imagesize(infile) maxwidth = width - 1 maxheight = height - 1 @@ -523,7 +524,7 @@ elif zoom <= -1: #}}} dfov = numpy.divide( - numpy.multpily(180, numpy.arctan(36/24)), numpy.pi) + numpy.multiply(180, numpy.arctan(36/24)), numpy.pi) pfact = 1 if pef == 0: # FIXME pfact = numpy.divide(0.01, dfov) @@ -541,7 +542,11 @@ 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)) + numpy.add( + numpy.multiply(width, width), + numpy.multiply(height, height) + ) + ) focal = numpy.divide(diag, numpy.multiply(2, tfov)) # calculate forward transform matrix Q @@ -636,6 +641,10 @@ for index in range(0, 4): vmax = vArr[index] delu = numpy.add(numpy.subtract(umax, umin), 1) delv = numpy.add(numpy.subtract(vmax, vmin), 1) +#FIXME ? +offsetu = 0 +offsetv = 0 +#added above if auto == "c": offsetu = numpy.divide(numpy.subtract(width, delu), 2) offsetv = numpy.divide(numpy.subtract(height, delv), 2) @@ -654,37 +663,37 @@ elif auto == "zc": offsetu = 0 offsetv = numpy.divide( numpy.subtract( - height, numpy.divide(numpy.multpily(delv, width), delu)), 2) + height, numpy.divide(numpy.multiply(delv, width), delu)), 2) else: _del = delv offsetu = numpy.divide( numpy.subtract( - width, numpy.divide(numpy.multpily(delu, height), delv)), 2) + width, numpy.divide(numpy.multiply(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))) + u1 = numpy.add( + offsetu, + numpy.multiply(numpy.subtract(u1, umin), numpy.divide(width, _del))) + v1 = numpy.add( + offsetv, + numpy.multiply(numpy.subtract(v1, vmin), numpy.divide(height, _del))) + u2 = numpy.add( + offsetu, + numpy.multiply(numpy.subtract(u2, umin), numpy.divide(width, _del))) + v2 = numpy.add( + offsetv, + numpy.multiply(numpy.subtract(v2, vmin), numpy.divide(height, _del))) + u3 = numpy.add( + offsetu, + numpy.multiply(numpy.subtract(u3, umin), numpy.divide(width, _del))) + v3 = numpy.add( + offsetv, + numpy.multiply(numpy.subtract(v3, vmin), numpy.divide(height, _del))) + u4 = numpy.add( + offsetu, + numpy.multiply(numpy.subtract(u4, umin), numpy.divide(width, _del))) + v4 = numpy.add( + offsetv, + numpy.multiply(numpy.subtract(v4, vmin), numpy.divide(height, _del))) # # now do the perspective distort if auto == "out": -- cgit v1.2.3-70-g09d2