#!/usr/bin/perl

use warnings;
use strict;

# This script takes a scan of a Kodak Disk and split it into 15
# cropped and properly rotated images.
#
# 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 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

######################################################################
# Converted 2015/7/2 to Perl v5 by Joe Touch / touch@isi.edu 
# Code changed from original ksh with the help of sh2p.pl
# new comments prefixed by "Touch:"
# description above edited to remove references to ksh
######################################################################

# 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

# Touch: here is a list of all output types and resulting arguments

my %outargs = ("jpg", "",
	       "tif", "-compress LZW -verbose");

# Touch: changed default to JPEG
my $Type = "jpg"; # or "tif"

my $Extra = $outargs{$Type};

# You don't have to change anything below here
my $Debug = 0;
# Normaly set to 0 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 !
my $RadiusMax = 8000;
my $DestDir = "";

my $argcount = scalar(@ARGV);

if (($argcount != 7) && ($argcount != 8)) {
    my $diemsg = "Usage: $0 X1 Y1 X6 Y6 X11 Y11 ImageFile [DestDir]\n" . 
	" Where Xn Yn are the coordinates of the point \n" . 
	" where the two images near '^N' do touch each others.\n" . 
	" If DestDir is not given, the images are saved in the current directory.\n" .
	"\n" .
	"and you gave $argcount arguments.\n";
    die $diemsg;
}

if ($Debug) {
    print "Argument list is:\n";
    print "Xa = $ARGV[0]\n";
    print "Ya = $ARGV[1]\n";
    print "Xb = $ARGV[2]\n";
    print "Yb = $ARGV[3]\n";
    print "Xc = $ARGV[4]\n";
    print "Yc = $ARGV[5]\n";
    print "Source = $ARGV[6]\n";
    if ($argcount == 8) {
	print "DestDir = $ARGV[7]\n\n";
    } else {
	print "DestDir = \n\n";
    }
}

# Touch: check output directory and create it if missing
if ( $argcount == 8 ) {
    $DestDir = $ARGV[7];
    if (!(-d $DestDir)) {
	mkdir($DestDir, 0644) || die "Cannot create directory $DestDir";
    }
    if (!(-w $DestDir)) {
	die "Cannot write to directory $DestDir";
    }
}

#  Given three points in space (A,B,C)
my $Xa = $ARGV[0];
my $Ya = $ARGV[1];
 # Point 1
my $Xb = $ARGV[2];
my $Yb = $ARGV[3];
# Point 6
my $Xc = $ARGV[4];
my $Yc = $ARGV[5];
# Point 11

my @argname = ('X1','Y1','X6','Y6','X11','Y11');

# Touch: check that the inputs are numeric
for (my $i = 0; $i < 6; $i++) {
    if (!($ARGV[$i] =~ /^\d+$/)) {
	die "Argument $argname[$i] = $ARGV[$i] must be numeric.";
    }
}

my $Source = $ARGV[6];
# Touch: input file

# Touch: check that it input file exists and is accessible
if (!(-e $Source)) {
    die "Input file $Source doesn't exist.";
} elsif (!(-r $Source)) {
    die "Input file $Source exists but isn't readable.";
}

my $FileName = "";

if ($Source =~ /([^\/]+)\.\w+$/) {
    $FileName = $1;
} else {
    die "Source file $Source missing a file suffix.\n";
}

my $Ratio= 1.30;
# Image ratio
# $R2D= 180./acos(-1);
# 180/pi, * for radian to degree conversion, / for the opposite
# Touch: use the Perl package instead, which includes rad2deg():

use Math::Trig;

# 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
my $AB = sqrt( ($Xa - $Xb) * ($Xa - $Xb) + ($Ya - $Yb) * ($Ya - $Yb) );
my $BC = sqrt( ($Xb - $Xc) * ($Xb - $Xc) + ($Yb - $Yc) * ($Yb - $Yc) );
my $AC = sqrt( ($Xa - $Xc) * ($Xa - $Xc) + ($Ya - $Yc) * ($Ya - $Yc) );

if ( $Debug ) {
    print "AB = $AB\n";
    print "BC = $BC\n";
    print "AC = $AC\n";
}

#   Direction cosines of AB(ABi,ABj) and AC(ACi,ACj)
my $ABi = ($Xb - $Xa) / $AB;
my $ABj = ($Yb - $Ya) / $AB;
my $ACi = ($Xc - $Xa) / $AC;
my $ACj = ($Yc - $Ya) / $AC;

if ( $Debug ) {
    print "Direction cosine ABi= $ABi\n";
    print "Direction cosine ABj= $ABj\n";
    print "Direction cosine ACi= $ACi\n";
    print "Direction cosine ACj= $ACj\n";
}

#   Cosine of angle BAC
my $cosBAC = ($AB * $AB + $AC * $AC - $BC * $BC) / (2.0 *  $AB * $AC);
my $AD = $cosBAC * $AC;
my $CD = sqrt($AC * $AC - $AD * $AD);
if ( $Debug ) {
    print "AD = $AD\n";
    print "CD = $CD\n";
}

#   Position of point D, which is C projected normally onto AB
my $Xd = $Xa + ($AD *  $ABi);
my $Yd = $Ya + ($AD *  $ABj);

if ( $Debug ) {
    print "Xd = $Xd\n";
    print "Yd = $Yd\n";
}

#   Direction cosines of CD(Cdi,CDj)
my $CDi = ($Xc - $Xd) / $CD;
my $CDj = ($Yc - $Yd) / $CD;

if ( $Debug ) {
    print "CDi = $CDi\n";
    print "CDj = $CDj\n";
}

# 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
my $sinBAC = sqrt(1.0 - $cosBAC * $cosBAC);
my $R = $BC / ($sinBAC * 2.0);

if ( $Debug ) {
    print "sinBAC = $sinBAC\n";
}

#   Centre of circumscribed circle is point E
my $X2e = $AB / 2.0;
my $Y2e = sqrt($R * $R - $X2e * $X2e);

if ( $Debug ) {
    print "X2e = $X2e\n";
    print "Y2e = $Y2e\n";
}

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

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

# Rotation to apply to get 1st image aligned properly on the right:
my $BaseRot = atan2($Ya - $Y, $Xa - $X ) - deg2rad(12.0);

if ( $Debug ) {
    print "Center X=$X\n";
    print "Center Y=$Y\n";
    print "BaseRot = ", rad2deg($BaseRot), "deg\n";
}

# Position of image to crop (after rotation)
my $Height = 2.0 * $R * sin(deg2rad(12.0));
my $Width = int($Height * $Ratio);
$Height= int($Height);

# Half the diagonal of the image
my $ImgRad = sqrt($Width * $Width + $Height * $Height) / 2.0;
my $r = $R * cos(deg2rad(12.0));
# Distance between circle center and point P (middle of inner 1st image height)
my $Rimg = $r + $Width / 2.0;
# Distance between circle center and middle of the images
my $IR = int($ImgRad);
my $IR2 = int($ImgRad * 2.0);
print "ImgRad = $ImgRad < Height = $Height < Width = $Width < r = ", int($r), 
    " < R = ", int($R), " < Rimg = ", int($Rimg), "\n";

if ( int($Rimg) > $RadiusMax ) {
    die "Resulting radius too high. Check your coordinates or increase RadiusMax\n";
}

# Prepare command
my $Exe = "convert $Source ";

# 15 rotations of 360/15=24 degrees around center
for (my $i = 1; $i <= 15; $i++) {
    # Rotation of the image
    my $Rot = $BaseRot + ($i - 1) * deg2rad(24.0);
    # Center of the image
    my $Ximg = $X + $Rimg * cos($Rot);
    my $Yimg = $Y + $Rimg * sin($Rot);
    # Adjust for out of frame crop
    my $FirstCrop = $IR2 . "x" . $IR2 . "+" . int($Ximg - $ImgRad) . "+" . int($Yimg - $ImgRad);
    # Touch: change any +- to -
    $FirstCrop =~ s/\+\-/\-/g;
    my $SRT = "$IR,$IR\\ " . (360 - rad2deg($Rot));
    my $DestFile = $DestDir . $FileName . "-KD" . $i . "\.$Type";
    if ( $Debug ) {
        print "$DestFile , $FirstCrop , $SRT\n";
    }

# 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  or the whole thing

    my $oneshotargs = 
	"-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 ";

#    my $doitnow = $Exe . $oneshotargs . "$DestFile";
#    print "NOW = $doitnow\n\n";
#    system($doitnow);
#    die "stopped after the first one";

    if ( $i < 15 ) {
	$Exe .= 
	    "\\( " . 
	    "+clone " .
	    $oneshotargs . "-write $DestFile " .
	    "+delete " .
	    "\\) ";
    } else {
	$Exe .= $oneshotargs . "$DestFile";
    }

}

if ( $Debug ) {
    print "Exe = $Exe\n";
} 

# Execute the whole thing
system "$Exe";