#! /bin/ksh
# This script takes a scan of a Kodak Disk and split it into 15 cropped and properly rotated images.
# You need the Korn shell (ksh, pdksh or mksh) for floating points arithmetics
#
# First open your image with a graphic editor and find the X/Y coordinates of the following 3 points:
# - The point (1) where the white parts of the images 1 and 2 touch each other (in front of ^1)
# - The point (6) in front of ^6 (or any other point between images)
# - The point (11) in front of ^11 (or any other point between images)
# You can use a partial scan and pick any points where images touch
# 
# Source: the maths were simplified from http://www.eng-tips.com/viewthread.cfm?qid=248118&page=2
# Caveat: I haven't tested it in Windows, but it should work in Cygwin as long as you have ksh and ImageMagick
#
# (c) 2009/11/06 Guillaume Dargaud, free use and modification by anyone
#     2009/11/16 Optimized for speed
# Tutorial: http://www.gdargaud.net/Hack/KodakDiskScan.html

# Those are the only things you (normally) should change in this script
# try "identify -list format | grep rw" to know which file formats you can write to
# and the optional -compress option for the corresponding file format
# I then post-process the resulting TIF files in a graphic editor
Type="tif"
Extra="-compress LZW -verbose"

# Other possibility:
#Type="jpg"
#Extra=""


# You don't have to change anything below here
Debug=""		# Normaly set to "" or 1

# Will stop the script if the circle diameter goes above this number of pixels
# In case you mistyped the coordinates, it will avoid having the script
# request 10Gb for a 160000pix image and the PC thrashing all night !
RadiusMax=8000

if [ $# -eq 8 ]; 
then
	DestDir=$8/
else
	if [ $# -eq 7 ]; 
	then
		DestDir=""
	else
		echo 1>&2 "Usage: $0 X1 Y1 X6 Y6 X11 Y11 ImageFile [DestDir]"
		echo 1>&2 "Where Xn Yn are the coordinates of the point where the two images near '^N' do touch each others."
		echo 1>&2 "If DestDir is not given, the images are saved in the current directory"
		exit 127
	fi
fi

#  Given three points in space (A,B,C)
Xa=$1;  Ya=$2	# Point 1
Xb=$3;  Yb=$4	# Point 6
Xc=$5;  Yc=$6	# Point 11
Source=$7
Suffix=$( echo $Source | sed -e "s/.*\.//" )
FileName=$( basename $Source .$Suffix )

(( Ratio=1.30 ))	# Image ratio
(( R2D=180./acos(-1) ))	# 180/pi, * for radian to degree conversion, / for the opposite

# Result of a formula in bash see "man bc"):
# A=$(echo "..." | bc -l)
# Result of a formula in ksh
# (( A=... ))

#   Lengths of AB, BC, AC
(( AB=sqrt( (Xa-Xb)*(Xa-Xb) + (Ya-Yb)*(Ya-Yb) ) ))
(( BC=sqrt( (Xb-Xc)*(Xb-Xc) + (Yb-Yc)*(Yb-Yc) ) ))
(( AC=sqrt( (Xa-Xc)*(Xa-Xc) + (Ya-Yc)*(Ya-Yc) ) ))
if [ $Debug ]; then echo "AB=$AB"; echo "BC=$BC"; echo "AC=$AC"; fi

#   Direction cosines of AB(ABi,ABj) and AC(ACi,ACj)
(( ABi=(Xb - Xa) / AB ))
(( ABj=(Yb - Ya) / AB ))
(( ACi=(Xc - Xa) / AC ))
(( ACj=(Yc - Ya) / AC ))
if [ $Debug ]; then echo "Direction cosine ABi=$ABi"; echo "Direction cosine ABj=$ABj"; echo "Direction cosine ACi=$ACi"; echo "Direction cosine ACj=$ACj"; fi

#   Cosine of angle BAC
(( cosBAC=(AB*AB + AC*AC - BC*BC) / (2. * AB * AC) ))
(( AD=cosBAC * AC ))
(( CD=sqrt(AC*AC - AD*AD) ))
if [ $Debug ]; then echo "AD=$AD"; echo "CD=$CD"; fi

#   Position of point D, which is C projected normally onto AB
(( Xd=Xa + (AD * ABi) ))
(( Yd=Ya + (AD * ABj) ))
if [ $Debug ]; then echo "Xd=$Xd"; echo "Yd=$Yd"; fi

#   Direction cosines of CD(Cdi,CDj)
(( CDi=(Xc - Xd) / CD ))
(( CDj=(Yc - Yd) / CD ))
if [ $Debug ]; then echo "CDi=$CDi"; echo "CDj=$CDj"; fi

# Diameter of circumscribed circle of a triangle is equal to
#   the length of any side divided by sine of the opposite angle.
#   This is done in a coordinate system where X is colinear with AB, Y is // to CD,
#   and Z is the normal (N) to X and Y, and the origin is point A
#         R = D / 2
(( sinBAC=sqrt(1. - cosBAC*cosBAC) ))
(( R= BC / (sinBAC*2) ))
if [ $Debug ]; then echo "sinBAC=$sinBAC"; fi

#   Centre of circumscribed circle is point E
(( X2e=AB / 2 ))
(( Y2e=sqrt(R*R - X2e*X2e) ))
if [ $Debug ]; then echo "X2e=$X2e"; echo "Y2e=$Y2e"; fi

#   Transform matrix
#                    Rotations                 Translations
#               ABi  ,   ABj  ,  ABk                 Xa
#               CDi  ,   CDj  ,  CDk                 Ya
#                Ni  ,    Nj  ,   Nk                 Za

#   Position of circle centre in absolute axis system
(( X=Xa + X2e*ABi + Y2e*CDi ))
(( Y=Ya + X2e*ABj + Y2e*CDj ))

# Rotation to apply to get 1st image aligned properly on the right:
(( BaseRot=atan2( Ya-Y, Xa-X )-12./R2D ))
if [ $Debug ]; then echo "Center X=$X"; echo "Center Y=$Y"; echo "BaseRot=$((R2D*BaseRot))deg"; fi

# Position of image to crop (after rotation)
(( Height=2*R*sin(12./R2D) ))
(( Width=int(Height*Ratio) ))
(( Height=int(Height) ))

# Half the diagonal of the image
(( ImgRad=sqrt(Width*Width+Height*Height)/2 ))
(( r=R*cos(12./R2D) ))		# Distance between circle center and point P (middle of inner 1st image height)
(( Rimg=r+Width/2 ))		# Distance between circle center and middle of the images
(( IR=int(ImgRad) ))
(( IR2=int(ImgRad*2) ))
echo "ImgRad=$ImgRad < Height=$Height < Width=$Width < r=$((int(r))) < R=$((int(R))) < Rimg=$((int(Rimg)))"

if [[ $((int($Rimg))) -gt $RadiusMax ]]
then
	echo "Resulting radius too high. Check your coordinates or increase RadiusMax"
	exit
fi

# Prepare command
Exe="convert $Source"

# 15 rotations of 360/15=24 degrees around center
for i in $( seq -w 1 15 )
do
	# Rotation of the image	
	(( Rot=BaseRot+(i-1)*24/R2D ))

	# Center of the image
	(( Ximg=X+Rimg*cos(Rot) ))
	(( Yimg=Y+Rimg*sin(Rot) ))
	# Adjust for out of frame crop
	FirstCrop=$( echo ${IR2}x${IR2}+$((int(Ximg-ImgRad)))+$((int(Yimg-ImgRad))) | sed -e "s/+-/-/g" )
	SRT="$IR,$IR $((360-R2D*Rot))"

	DestFile=${DestDir}$FileName-KD$i.$Type
	if [ $Debug ]; then echo "$DestFile, $FirstCrop, $SRT"; fi

	# This takes only 10~20s per image because it does a crop before the rotation and the final precise crop
#	nice convert $Source \
#		-crop $FirstCrop\!	\
#		-background lightblue -flatten +repage \
#		-distort SRT "$SRT"	\
#		-crop $((int(Width)))x$((int(Height)))+$((int(ImgRad-Width/2)))+$((int(ImgRad-Height/2)))\! \
#		+repage \
#		$Extra $DestFile

	# Same as above but we pack all 15 extractions into one single command: 40s for the whole thing
	if [[ $i -lt 15 ]]
	then
		Exe="$Exe \( +clone \
			-crop $FirstCrop\!	\
			-background lightblue -flatten +repage \
			-distort SRT \"$SRT\"	\
			-crop $((int(Width)))x$((int(Height)))+$((int(ImgRad-Width/2)))+$((int(ImgRad-Height/2)))\! \
			+repage \
			$Extra -write $DestFile +delete \) "
	else
		Exe="$Exe -crop $FirstCrop\!	\
			-background lightblue -flatten +repage \
			-distort SRT \"$SRT\"	\
			-crop $((int(Width)))x$((int(Height)))+$((int(ImgRad-Width/2)))+$((int(ImgRad-Height/2)))\! \
			+repage \
			$Extra $DestFile"
	fi
done

if [ $Debug ]; then echo "$Exe"; fi

# Execute the whole thing
eval "$Exe"