#!/bin/bash
# 
# Developed by Fred Weinhaus 4/14/2008 .......... revised 4/25/2015
#
# ------------------------------------------------------------------------------
# 
# Licensing:
# 
# Copyright © Fred Weinhaus
# 
# My scripts are available free of charge for non-commercial use, ONLY.
# 
# For use of my scripts in commercial (for-profit) environments or 
# non-free applications, please contact me (Fred Weinhaus) for 
# licensing arrangements. My email address is fmw at alink dot net.
# 
# If you: 1) redistribute, 2) incorporate any of these scripts into other 
# free applications or 3) reprogram them in another scripting language, 
# then you must contact me for permission, especially if the result might 
# be used in a commercial or for-profit environment.
# 
# My scripts are also subject, in a subordinate manner, to the ImageMagick 
# license, which can be found at: http://www.imagemagick.org/script/license.php
# 
# ------------------------------------------------------------------------------
# 
####
#
# USAGE: defisheye [-i ifov] [-o ofov] [-t type] [-f format] [-c xc,yc] [-r radius] [-a angle] infile outfile
# USAGE: defisheye -d [-c xc,yc] [-r radius] [-s strokecolor] infile outfile
# USAGE: defisheye [-h or -help]
# 
# OPTIONS:
# 
# -i      ifov           input (fisheye) image field of view in degrees;
#                        float; 0<ifov<=180; default=180
# -o      ofov           output (perspective) image field of view in degrees;
#                        float; 0<ifov<180; default=120
# -t      type           type of fisheye lens; linear, equalarea, 
#                        orthographic, stereographic; default=linear
# -f      format         format of fisheye lens image; circular or fullframe;
#                        default=circular
# -c      xc,yc          coordinates of the center of the fisheye image area, 
#                        which will also become the center of the perspective 
#                        view; float; default=center of fisheye image
# -r      radius         radius of fisheye area in the input image; float;
#                        radius>0; default=min(width,height)/2
# -a      angle          angle of rotation (clockwise) for the perspective image;
#                        float; 0<=angle<360; default=0
# -d                     draw circle on the input fisheye image of specified radius, 
#                        center and color=scolor, which will then become the output
# -s      scolor         scolor is the stroke color for a circle that will be 
#                        drawn on the input fisheye image; default=white
# 
###
# 
# NAME: DEFISHEYE 
# 
# PURPOSE: To correct for fisheye distortion in an image.
# 
# DESCRIPTION: DEFISHEYE is designed to transform a fisheye image into 
# a normal perspective view looking towards the center of the fisheye image.
# 
# 
# ARGUMENTS: 
# 
# -i ifov ... IFOV is the input fisheye image field of view in degrees. A 
# value of 180 will correspond to a hemispherical fisheye image within the 
# circular area. Values are floats in the range 0<ifov<=180. The default is 
# 180 degrees for a full hemisphere.
# 
# -o ofov ... OFOV is the output perspective image field of view in
# degrees. Values are floats in the range 0<ofov<180. The default is to
# use 120 degrees both vertically and horizontally for a circular fisheye
# and diagonally for a full frame fisheye. The value for ofov relative to
# the ifov determines the proportional amount of the fisheye area that
# will be transformed. Note that in comparison, a value of 48.8 degrees
# corresponds to a diagonal field of view from a 35 mm camera (film size
# 36mm x 24mm) with a 50mm focal length lens, i.e. a "normal" view.
# Similarly, when the image diagonal is equal to the focal length of the
# camera, the field of view is about 53.1 degrees. Although the default
# value is perhaps not appropriate to a normal perspective image, this
# will produce an image that maximizes the area covered, but without too
# much distortion. If the original fisheye image was viewed obliquely,
# i.e. the camera was tilted between horizontal and vertical, then the
# resulting perspective view will have perspective distortion. That is
# vertical edges will be tilted outward or inward. Post processing with my
# 3Drotate or rotate3D script will then correct for this perspective
# distortion. Perspective distortion will be more pronounced with larger
# values for ofov.
# 
# -t type ... TYPE is the type of fisheye lens. The choices are: linear 
# (equidistant), equalarea (equisolid), orthographic and stereographic. 
# The default is linear.
# 
# -f format ... FORMAT is the format of the fisheye lens image. The choices are: 
# circular (image fills a circle that spans the minimum of the width or height) 
# or fullframe (image spans a circle that spans the diagonal dimension). The 
# default is circular.
# 
# -c xc,yc ... XC,YC are the pixel coordinates in the input fisheye  
# image that correspond to the (circular fisheye area) center. The pixel at 
# this coordinate will then become the center of the perspective image. The
# default values are the center of the input fisheye image. Values  
# are non-negative floats. You can use the -d option to validate your 
# choice of center and radius for the fisheye image. See more below.
# 
# -r radius ... RADIUS is the radius of the fisheye circular area in the 
# input image. Values are floats greater than zero. The default is half the 
# minimum value between the input image width and height.
# 
# -a angle ... ANGLE is the clockwise positive rotation angle for the output 
# perspective image relative to the orientation of the input fisheye image.  
# Values are non-negative floats in range 0<=angle<360. The default is 0.
# 
# -d ... Use of this argument produces an ouput image that is simply the input 
# image with a circle drawn on it to show where the expected fisheye image 
# area is located. You can specify a radius and center point if you want to 
# adjust the transformation to use the precise center and radius that matches 
# the area delimited by the circle. Radius and center default as described 
# above.
# 
# -s scolor ... SCOLOR is the stroke color to use to draw the circle when the 
# -d argument is used. The default is white.
# 
# See the following references for definitions and mathematical details of each 
# type of fisheye lens: 
# http://en.wikipedia.org/wiki/Fisheye_lens
# http://www.bobatkins.com/photography/technical/field_of_view.html
# 
# 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. 
# 
######
# 

# set default values
ifov=180					# fisheye field of view (aperture) in degrees
ofov=120					# perspective field of view (aperture) in degrees
xc=""						# center of fisheye area
yc=""						# center of fisheye area
rad=""						# radius of fisheye area
ang=0						# image rotation in degrees clockwise
type="linear"				# linear, equalarea, orthographic, stereographic
format="circular"			# circular, fullframe
draw="no"					# flag to draw circle on input
sc="red"					# stroke color for drawing circle

# set directory for temporary files
dir="."    # suggestions are dir="." or dir="/tmp"

# set up functions to report Usage and Usage with Description
PROGNAME=`type $0 | awk '{print $3}'`  # search for executable on path
PROGDIR=`dirname $PROGNAME`            # extract directory of program
PROGNAME=`basename $PROGNAME`          # base name of program
usage1() 
	{
	echo >&2 ""
	echo >&2 "$PROGNAME:" "$@"
	sed >&2 -e '1,/^####/d;  /^###/g;  /^#/!q;  s/^#//;  s/^ //;  4,$p' "$PROGDIR/$PROGNAME"
	}
usage2() 
	{
	echo >&2 ""
	echo >&2 "$PROGNAME:" "$@"
	sed >&2 -e '1,/^####/d;  /^######/g;  /^#/!q;  s/^#*//;  s/^ //;  4,$p' "$PROGDIR/$PROGNAME"
	}


# function to report error messages
errMsg()
	{
	echo ""
	echo $1
	echo ""
	usage1
	exit 1
	}


# function to test for minus at start of value of second part of option 1 or 2
checkMinus()
	{
	test=`echo "$1" | grep -c '^-.*$'`   # returns 1 if match; 0 otherwise
    [ $test -eq 1 ] && errMsg "$errorMsg"
	}

# test for correct number of arguments and get values
if [ $# -eq 0 ]
	then
	# help information
   echo ""
   usage2
   exit 0
elif [ $# -gt 16 ]
	then
	errMsg "--- TOO MANY ARGUMENTS WERE PROVIDED ---"
else
	while [ $# -gt 0 ]
		do
			# get parameter values
			case "$1" in
		  -h|-help)    # help information
					   echo ""
					   usage2
					   exit 0
					   ;;
				-i)    # get ifov
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID IFOV SPECIFICATION ---"
					   checkMinus "$1"
					   ifov=`expr "$1" : '\([.0-9]*\)'`
					   [ "$ifov" = "" ] && errMsg "--- IFOV=$ifov MUST BE A NON-NEGATIVE FLOAT ---"
					   ifovtestA=`echo "$ifov <= 0" | bc`
					   ifovtestB=`echo "$ifov > 180" | bc`
					   [ $ifovtestA -eq 1 -o $ifovtestB -eq 1 ] && errMsg "--- IFOV=$ifov MUST BE A FLOAT GREATER THAN 0 AND LESS THAN OR EQUAL 180 ---"
					   ;;
				-o)    # get ofov
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID OFOV SPECIFICATION ---"
					   checkMinus "$1"
					   ofov=`expr "$1" : '\([.0-9]*\)'`
					   [ "$ofov" = "" ] && errMsg "--- IFOV=$ofov MUST BE A NON-NEGATIVE FLOAT ---"
					   ofovtestA=`echo "$ofov <= 0" | bc`
					   ofovtestB=`echo "$ofov >= 180" | bc`
					   [ $ofovtestA -eq 1 -o $ofovtestB -eq 1 ] && errMsg "--- OFOV=$ofov MUST BE A FLOAT GREATER THAN 0 AND LESS THAN 180 ---"
					   ;;
				-t)    # get  type
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID TYPE SPECIFICATION ---"
					   checkMinus "$1"
					   type="$1"
					   [ "$type" != "linear" -a "$type" != "equalarea" -a "$type" != "orthographic" -a "$type" != "stereographic" ] && errMsg "--- INVALID TYPE VALUE ---"
					   ;;
				-f)    # get  format
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID FORMAT SPECIFICATION ---"
					   checkMinus "$1"
					   format="$1"
					   [ "$format" != "circular" -a "$format" != "fullframe" ] && errMsg "--- INVALID TYPE VALUE ---"
					   ;;
				-c)    # get coords
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID COORDS SPECIFICATION ---"
					   checkMinus "$1"
					   test=`echo "$1" | tr "," " " | wc -w`
					   [ $test -eq 1 -o $test -gt 2 ] && errMsg "--- INCORRECT NUMBER OF COORDINATES SUPPLIED ---"
					   coords=`expr "$1" : '\([.0-9]*,[.0-9]*\)'`
					   [ "$coords" = "" ] && errMsg "--- COORDS=$coords MUST BE A PAIR OF NON-NEGATIVE FLOATS SEPARATED BY A COMMA ---"
					   coords="$1,"
		   			   xc=`echo "$coords" | cut -d, -f1`
		   			   yc=`echo "$coords" | cut -d, -f2`
					   ;;
				-r)    # get rad
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID RADIUS SPECIFICATION ---"
					   checkMinus "$1"
					   rad=`expr "$1" : '\([.0-9]*\)'`
					   [ "$rad" = "" ] && errMsg "--- RADIUS=$rad MUST BE A NON-NEGATIVE FLOAT ---"
					   radtestA=`echo "$rad <= 0" | bc`
					   [ $radtestA -eq 1 ] && errMsg "--- RADIUS=$rad MUST BE A FLOAT GREATER THAN 0 ---"
					   ;;
				-a)    # get ang
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID ANGLE SPECIFICATION ---"
					   checkMinus "$1"
					   ang=`expr "$1" : '\([.0-9]*\)'`
					   [ "$ang" = "" ] && errMsg "--- ANGLE=$ang MUST BE A NON-NEGATIVE FLOAT ---"
					   angtestA=`echo "$ang < 0" | bc`
					   angtestB=`echo "$ang >= 360" | bc`
					   [ $angtestA -eq 1 -o $angtestB -eq 1 ] && errMsg "--- ANGLE=$ang MUST BE A FLOAT GREATER THAN OR EQUAL TO 0 AND LESS THAN 360 ---"
					   ;;
				-s)    # get  scolor
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID STROKE COLOR SPECIFICATION ---"
					   checkMinus "$1"
					   sc="$1"
					   ;;
				-d)    # enable draw circle
					   draw="yes"
					   ;;
				 -)    # STDIN and end of arguments
					   break
					   ;;
				-*)    # any other - argument
					   errMsg "--- UNKNOWN OPTION ---"
					   ;;
		     	 *)    # end of arguments
					   break
					   ;;
			esac
			shift   # next option
	done
	#
	# get infile and outfile
	infile="$1"
	outfile="$2"
fi

# test that infile provided
[ "$infile" = "" ] && errMsg "NO INPUT FILE SPECIFIED"

# test that outfile provided
[ "$outfile" = "" ] && errMsg "NO OUTPUT FILE SPECIFIED"

tmpA="$dir/defisheye_$$.mpc"
tmpB="$dir/defisheye_$$.cache"
trap "rm -f $tmpA $tmpB;" 0
trap "rm -f $tmpA $tmpB; exit 1" 1 2 3 15
trap "rm -f $tmpA $tmpB; exit 1" ERR

# get im_version
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" -ge "07000000" ]; then
	identifying="magick identify"
else
	identifying="identify"
fi

if convert -quiet "$infile" +repage "$tmpA"
	then
	width=`$identifying -format %w $tmpA`
	height=`$identifying -format %h $tmpA`
	
	# compute dimension depending upon format
	if [ "$format" = "circular" ]
		then
		dim=`convert xc: -format "%[fx:min($width,$height)]" info:`
	elif [ "$format" = "fullframe" ]
		then
		dim=`convert xc: -format "%[fx:hypot($width,$height)]" info:`
	fi
	[ "$rad" != "" ] && dim=`convert xc: -format "%[fx:2*$rad]" info:`


	# compute half-widths
	w2=`convert xc: -format "%[fx:($width-1)/2]" info:`
	h2=`convert xc: -format "%[fx:($height-1)/2]" info:`
	
	# compute center of image
	if [ "$xc" = "" ]
		then
		xcd=$w2
		xcs=$xcd
	else
		xcd=$w2
		xcs=$xc
	fi
	if [ "$yc" = "" ]
		then
		ycd=$h2
		ycs=$ycd
	else
		ycd=$h2
		ycs=$yc
	fi

	# compute output (perspective) focal length and its inverse from ofov
	# phi=fov/2; r=N/2
	# r/f=tan(phi); f=r/tan(phi); f= (N/2)/tan((fov/2)*(pi/180)) = N/(2*tan(fov*pi/360))
	ofoc=`convert xc: -format "%[fx:$dim/(2*tan($ofov*pi/360))]" info:`
	ofocinv=`convert xc: -format "%[fx:1/$ofoc]" info:`
	
	# compute sin(ang) and cos(ang)
	[ "$ang" != "0" -a "$ang" != "0.0" ] && cosang=`convert xc: -format "%[fx:cos($ang*pi/180)]" info:`
	[ "$ang" != "0" -a "$ang" != "0.0" ] && sinang=`convert xc: -format "%[fx:sin($ang*pi/180)]" info:`
	else
		errMsg "--- FILE $infile DOES NOT EXIST OR IS NOT AN ORDINARY FILE, NOT READABLE OR HAS ZERO SIZE ---"
fi

# Pertinent equations:
# note phi=fov/2; fov=field of view (aperture)
# note r=N/2; N=min(width,height)
# perspective: r=f*tan(phi); f=r/tan(phi); f=(N/2)/tan((fov/2)*(pi/180))=N/(2*tan(fov*pi/360))
# linear: r=f*phi; f=r/phi; f=(N/2)/((fov/2)*(pi/180))=N*180/(fov*pi)
# equalarea: r=2*f*sin(phi/2); f=(r/2)/sin(phi/2); f=(N/4)/(sin((fov/4)*(pi/180)))=N/(4*sin(fov*pi/720))
# orthographic: r=f*sin(phi); f=r/sin(phi); f=(N/2)/sin((fov/2)*(pi/180))=N/(2*sin(fov*pi/360))
# stereographic: r=2*f*tan(phi/2); f=(r/2)/tan(phi/2); f=(N/4)/(tan((fov/4)*(pi/180)))=N/(4*tan(fov*pi/720))

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`

# see Bourke diagram at http://local.wasp.uwa.edu.au/~pbourke/projection/fish2/
# assume looking straight down on hemisphere
# convert i,j in output to rd,theta (where theta is same angle in input and output as output is looking straight down)
# then convert rd to phi using perspective equations
# then convert phi to rr (radius in input) using fisheye equations
# then convert rr,theta to x,y in input (note shortcut that avoids theta)

# define transform terms

# get rd (and theta unused) in output (destination) from i,j
xd="xd=i-$xcd;"
yd="yd=j-$ycd;"
if [ "$im_version" -ge "06030600" ]
	then 
	rd="rd=hypot(xd,yd);"
else
	rd="rd=sqrt(xd^2+yd^2);"
fi
: ' unused - see notes below
if [ "$ang" != "0" -a "$ang" != "0.0" ]
	then
	theta="theta=atan2(yd,xd)-$ang;"
else
	theta="theta=atan2(yd,xd);"
fi
'

#get phi from perspective radius
# phi is now a reserved word in fx since Oct, 2011
# so change phi to phiang
phiang="phiang=atan($ofocinv*rd);"

#convert phi to rr (input radius) according to type of fisheye
if [ "$type" = "linear" ]
	then
	ifoc=`convert xc: -format "%[fx:($dim*180)/($ifov*pi)]" info:`
	rr="rr=$ifoc*phiang;"
elif [ "$type" = "equalarea" ]
	then
	#note ifoc2 rather than ifoc to save computations
	ifoc2=`convert xc: -format "%[fx:$dim/(2*sin($ifov*pi/720))]" info:`
	rr="rr=$ifoc2*sin(phiang/2);"
elif [ "$type" = "orthographic" ]
	then
	ifoc=`convert xc: -format "%[fx:$dim/(2*sin($ifov*pi/360))]" info:`
	rr="rr=$ifoc*sin(phiang);"
elif [ "$type" = "stereographic" ]
	then
	#note ifoc2 rather than ifoc to save computations
	ifoc2=`convert xc: -format "%[fx:$dim/(2*tan($ifov*pi/720))]" info:`
	rr="rr=$ifoc2*tan(phiang/2);"
fi

# convert rr (and theta unused due to shortcut) to x,y in input (source)
# allow optional final image rotation
if [ "$ang" != "0" -a "$ang" != "0.0" ]
	then
	# use slower method when want rotation
	# this is slow:  2m32s (128x128)
#	xs="xs=rr*cos(theta)+$xcs;"
#	ys="ys=rr*sin(theta)+$ycs;"
	# this is a bit faster and avoids theta:  1m23s  (128x128)
	xs="xs=(rd?rr/rd:0)*xd*$cosang + (rd?rr/rd:0)*yd*$sinang + $xcs;"
	ys="ys=-(rd?rr/rd:0)*xd*$sinang + (rd?rr/rd:0)*yd*$cosang + $ycs;"
else
	# more efficient (by about 2x) to do this as avoids theta (but prevents rotation)
	# cos(theta)=xd/rd; sin(theta)=yd/rd
	# this is faster: 1m12s (128x128)
	xs="xs=(rd?rr/rd:0)*xd+$xcs;"
	ys="ys=(rd?rr/rd:0)*yd+$ycs;"
fi

# process image
if [ "$draw" = "yes" ]
	then
	[ "$rad" = "" ] && rad=`convert xc: -format "%[fx:min($width,$height)/2]" info:`
	xr=`convert xc: -format "%[fx:$xcs+$rad]" info:`
	convert $tmpA -fill none -stroke "$sc" -draw "circle $xcs,$ycs $xr,$ycs" "$outfile"
else
	convert $tmpA -virtual-pixel black -monitor -fx \
		"$xd $yd $rd $phiang $rr $xs $ys u.p{xs,ys}" \
		"$outfile"
fi
exit 0