#!/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, 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) 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 # 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