#!/usr/bin/python3

# Converted to Python3 by Joe Touch / touch@isi.edu
# April 4, 2016

# 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

import sys, os, argparse, math


# Touch: here is a list of all output types and resulting arguments
outargs = {"jpg": "", "tif": "-compress LZW -verbose"}
    
# Touch: changed default to JPEG
Type = "jpg"
# or "tif"
    
# Touch: this appears to force a check that the Type is one of the outargs
try:
    Extra = outargs[Type]
except:
    sys.exit("invalid output type " + Type)
    
# You don't have to change anything below here
Debug = False
# Normaly set to False or True

# 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;

parser = argparse.ArgumentParser(usage='''\
Usage: split.py X1 Y1 X6 Y6 X11 Y11 ImageFile [DestDir]
    Where Xn Yn are the coordinates of the point
    where the two images near ^ touch each other
    If DestDir is not given, the images are saved in the current directory''')
parser.add_argument("Xa", help="X1", type=int)
parser.add_argument("Ya", help="Y1", type=int)
parser.add_argument("Xb", help="X6", type=int)
parser.add_argument("Yb", help="Y6", type=int)
parser.add_argument("Xc", help="X11", type=int)
parser.add_argument("Yc", help="Y11", type=int)
parser.add_argument("Source", help="ImageFile")
parser.add_argument("DestDir", help="DestDir", nargs='?', default="")
args = parser.parse_args()

# Touch: get rid of extra \r's due to Windows-isms
args.Source = args.Source.rstrip('\r')
args.DestDir = args.DestDir.rstrip('\r')

if Debug:
    print("Argument list is:");
    print("Xa =", args.Xa)
    print("Ya =", args.Xb)
    print("Xb =", args.Xb)
    print("Yb =", args.Yb)
    print("Xc =", args.Xc)
    print("Yc =", args.Yc)
    print("Source =", args.Source)
    print("DestDir =", args.DestDir)
    print()

# Touch: check output directory and create it if missing
if (args.DestDir != ""):
    if not os.path.exists(args.DestDir):
        try:
            os.makedirs(args.DestDir)
        except:
            sys.exit("Cannot create directory " + args.DestDir)
    if not os.access(args.DestDir, os.W_OK):
        sys.exit("Cannot write to directory " + args.DestDir)

# Touch: check that it input file exists and is accessible
if not os.path.isfile(args.Source):
    sys.exit("Input file " + args.Source + " doesn't exist.")
elif not os.access(args.Source, os.R_OK):
    sys.exit("Input file " + args.Source + " exists but isn't readable.")

# Touch: extract filename prefix to use for each output image
FileName, FileExt = os.path.splitext(args.Source)

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

#   Lengths of AB, BC, AC
AB = math.sqrt( (args.Xa - args.Xb)**2 + (args.Ya - args.Yb)**2 )
BC = math.sqrt( (args.Xb - args.Xc)**2 + (args.Yb - args.Yc)**2 )
AC = math.sqrt( (args.Xa - args.Xc)**2 + (args.Ya - args.Yc)**2 )

if Debug:
    print("AB =", AB)
    print("BC =", BC)
    print("AC =", AC)
    
#   Direction cosines of AB(ABi,ABj) and AC(ACi,ACj)
ABi = (args.Xb - args.Xa) / AB;
ABj = (args.Yb - args.Ya) / AB;
ACi = (args.Xc - args.Xa) / AC;
ACj = (args.Yc - args.Ya) / AC;

if Debug:
    print("Direction cosine ABi=", ABi)
    print("Direction cosine ABj=", ABj)
    print("Direction cosine ACi=", ACi)
    print("Direction cosine ACj=", ACj)
    
#   Cosine of angle BAC
cosBAC = (AB**2 + AC**2 - BC**2) / (2.0 * AB * AC)
AD = cosBAC * AC
CD = math.sqrt(AC**2 - AD**2)

if Debug:
    print("AD =", AD)
    print("CD =", CD)
    
#   Position of point D, which is C projected normally onto AB
Xd = args.Xa + (AD * ABi)
Yd = args.Ya + (AD * ABj)

if Debug:
    print("Xd =", Xd)
    print("Yd =", Yd)

#   Direction cosines of CD(Cdi,CDj)
CDi = (args.Xc - Xd) / CD
CDj = (args.Yc - Yd) / CD

if Debug:
    print("CDi =", CDi)
    print("CDj =", CDj)

# 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 = math.sqrt(1.0 - cosBAC**2)
R = BC / (sinBAC * 2.0)

if Debug:
    print("sinBAC =", sinBAC)

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

if Debug:
    print("X2e =", X2e)
    print("Y2e =", Y2e)

#   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 = args.Xa + X2e * ABi + Y2e * CDi
Y = args.Ya + X2e * ABj + Y2e * CDj

# Rotation to apply to get 1st image aligned properly on the right:
BaseRot = math.atan2(args.Ya - Y, args.Xa - X ) - math.radians(12.0);

if Debug:
    print("Center X=", X)
    print("Center Y=", Y)
    print("BaseRot =", math.degrees(BaseRot),"deg")

# Position of image to crop (after rotation)
Height = 2.0 * R * math.sin(math.radians(12.0))
Width = math.floor(Height * Ratio)
Height= math.floor(Height)

# Half the diagonal of the image
ImgRad = math.sqrt(Width**2 + Height**2) / 2.0
r = R * math.cos(math.radians(12.0))
# Distance between circle center and point P (middle of inner 1st image height)
Rimg = r + Width / 2.0
# Distance between circle center and middle of the images
IR = math.floor(ImgRad)
IR2 = math.floor(ImgRad * 2.0)
print("ImgRad =", ImgRad, "< Height =", Height, "< Width =", Width, "< r =", math.floor(r), "< R =", math.floor(R), "< Rimg =", math.floor(Rimg))

if math.floor(Rimg) > RadiusMax:
    sys.exit("Resulting radius too high. Check your coordinates or increase RadiusMax")

# Prepare command
Exe = "convert " + args.Source + " "

# 15 rotations of 360/15=24 degrees around center
for i in range(1, 15+1):
    # Rotation of the image
    Rot = BaseRot + (i - 1) * math.radians(24.0)
    # Center of the image
    Ximg = X + Rimg * math.cos(Rot)
    Yimg = Y + Rimg * math.sin(Rot)
    # Adjust for out of frame crop
    FirstCrop = str(IR2) + "x" + str(IR2) + "+" + str(math.floor(Ximg - ImgRad)) + "+" + str(math.floor(Yimg - ImgRad))

    # Touch: change any +- to -
    FirstCrop.replace("\+\-","\-")
    SRT = str(IR) + "," + str(IR) + "\ " + str(360 - math.degrees(Rot))
    DestFile = args.DestDir.rstrip('/\r') + "/" + FileName + "-KD" + str(i) + "." + Type

    if Debug:
        print(DestFile, ",", FirstCrop, ",", SRT)

        # 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

        oneshotargs = "-crop " + FirstCrop + "! " + "-background lightblue -flatten +repage " + "-distort SRT " + SRT + " " + "-crop " + str(math.floor(Width)) + "x" + str(math.floor(Height)) + "+" + str(math.floor(ImgRad - Width / 2)) + "+" + str(math.floor(ImgRad - Height / 2)) + "! " + "+repage " + Extra + " " 
        
        #    doitnow = Exe + oneshotargs + DestFile;
        #    print("NOW = " + doitnow)
        #    os.system(doitnow);
        #    sys.exit("stopped after the first one")

        if i < 15:
            Exe += "\\( +clone " + oneshotargs + "-write " + DestFile + " +delete \\) "
        else:
            Exe += oneshotargs + DestFile

if Debug:
    print("Exe =", Exe)            

# Execute the whole thing
os.system(Exe)