path: root/bin/3Drotate.py

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