#! /usr/bin/env perl
use warnings;

#######################################################################
#
# Copyright(c) 2014 Whitehead Institute for Biomedical Research.
#              All Rights Reserved
#
#  Create a colored spiral based on genome/gene regions associated with colors
#
#  Written by Helen Skaletsky, Page Lab, Whitehead Institute
#  and George Bell, BaRC, Whitehead Institute
#
#  Version 1.0: 3 October 2014
#
############################################
#$ARGV[0] input file
#$ARGV[1] total length
#$ARGV[2] pixels per unit
#$ARGV[3] spiral width in pixels
#$ARGV[4] additional marks file (optional) => Same format as input file; will appear as lines next to spiral region
#
#######################################################################


if (! $ARGV[3])
{
	print STDERR "\nCreate a spiral postscript image from a vector of regions and colors\n";
	print STDERR "USAGE: $0 inputData regionLength resolution spiralWidth [marksData] > spiral.ps\n";
	print STDERR "Ex:	 $0 myRegions.txt 1000 6 8 > myRegions.spiral.ps\n";
	print STDERR "\nwhere inputData is a tab-delimited file containing columns for startPos, endPos, and color (as hex code).\n\n";
	exit;
}

$size = 600; 
$R = int(0.5 * $size);
$markWidth = 2;			# Width of line showing mark areas
$length = $ARGV[1]; 	# total length of the sequence in your units
$del = $ARGV[2]; 		# resolution: pixels/unit, recommended value is 6000/length
$wi = $ARGV[3]; 		# width of the spiral; 5-10 recommended
if ($ARGV[4])
{
	$marks = $ARGV[4];
}
$L = $length * $del;
$step = get_step($L, $R, $wi);

if($wi > 2*3.14*$step){
	print "Unprintable\n"; exit();
}

$title = '%!PS-Adobe-3.1';
print "$title\n";

open(IN, "$ARGV[0]") || die "Cannot open $ARGV[0] file\n";
<IN>;
$b = 0;
while(<IN>){
	chop; 
	($b1, $e1, $color) = split '\t';
	$color = fromHexToPScolor($color);
 
	$b1 *= $del; $e1 *= $del;
	if($b1 - $b > 0.5){$b1 = $b + 0.5;}
	$b = $b1; 
	$e = $b + 5; $f = 0;
	$d = 1;
	
	while($e < $e1 + 5){
		if($e > $e1){$e = $e1; $f = 1;}
		$s = coor($b, $step);
		($fi1, $r1) = split ' ', $s;
		$s = coor($e, $step);
		($fi2, $r2) = split ' ', $s;
		$x1 = $r1 * cos($fi1); $y1 = $r1 * sin($fi1);
		$x2 = $r2 * cos($fi2); $y2 = $r2 * sin($fi2);
		$r3 = $r1 + $wi; $r4 = $r2 + $wi;
		$x3 = $r3 * cos($fi1); $y3 = $r3 * sin($fi1);
		$x4 = $r4 * cos($fi2); $y4 = $r4 * sin($fi2);
		$x1 += $R; $x2 += $R; $y1 += $R; $y2 += $R;
		$x3 += $R; $x4 += $R; $y3 += $R; $y4 += $R;

		print_segment($x1, $y1, $x2, $y2, $x3, $y3, $x4, $y4, 0.1, $color ,$color);
		if($fi2 > 6){$d = 5;}
		$b = $e; $e += $d;
		if($f){last;}
	}

}
close IN;

if($marks){
	open(IN, $marks) || die "Cannot open marks file\n";
	while(<IN>){
		chop; 
		($b, $e, $color) = split '\t';
		$b--;	# GB
		$color = fromHexToPScolor($color);
		$b *= $del; $e *= $del; 
		$s = coor($b, $step);
		($fi1, $r1) = split ' ', $s; $r1 += $wi * 1.5;
		$s = coor($e, $step);
		($fi2, $r2) = split ' ', $s; 
		$fi = $fi1 + 0.01; $r = $step * $fi + $wi * 1.5;
		while($fi <= $fi2){
			$x1 = $r1 * cos($fi1); $y1 = $r1 * sin($fi1);
			$x = $r * cos($fi); $y = $r * sin($fi);
			$x1 += $R; $y1 += $R; $x += $R; $y += $R;
			$y1 += 180; $y += 180;
			print "newpath\n$markWidth setlinewidth\n";
			print "$color setrgbcolor\n$x $y moveto\n";
			print "$x1 $y1 lineto\nstroke\n";
			$fi1 = $fi; $fi += 0.01;
			$r1 = $step * $fi1 + 1.5 * $wi; 
			$r = $step * $fi + 1.5 * $wi;
		}
	}
	close IN;	
}

print "%EOF\n";	

sub print_segment{
	my($x1, $y1, $x2, $y2, $x3, $y3, $x4, $y4, $th, $gray, $color) = @_;
	print "newpath\n";

	$y1 += 180; $y2 += 180; $y3 += 180; $y4 += 180;

	print_line($x1, $y1, $x3, $y3, $th, $gray);
	print_line($x3, $y3, $x4, $y4, $th, $gray);
	print_line($x4, $y4, $x2, $y2, $th, $gray);
	print_line($x2, $y2, $x1, $y1, $th, $gray);
	print "closepath\ngsave\n";
	print "$color setrgbcolor fill\n";
	print "grestore\n";
	print "%EOF\n";
	
}

sub print_line{
	my($x1, $y1, $x2, $y2, $w, $color) = @_;

	print "$color setrgbcolor\n";
	print "$w setlinewidth\n";
	print "$x1 $y1 moveto\n$x2 $y2 lineto\n\n";
}

sub coor{
	my($dist, $a) = @_;
	my($fi0, $fi1, $del, $x0, $x1, $x2, $x3, $fi, $r, $t, $s);

	$fi0 = 0; $fi1 = 40*3.1415;

	while($fi1 - $fi0 > 0.001){
		$del = $fi1 - $fi0; $x0 = $fi0; $x1 = $x0 + $del/3; 
		$x2 = $x1 + $del/3; $x3 = $fi1;
			# $y0 = dist($x0, $a);  	# "main::y0" used only once [GB]
		$y1 = dist($x1, $a);
		$y2 = dist($x2, $a); $y3 = dist($x3, $a);
		if($dist < $y1){
			$fi1 = $x1; 
		}
		elsif($dist < $y2){
			$fi1 = $x2; $fi0 = $x1;
		}
		else{
			$fi1 = $x3; $fi0 = $x2;
		}
	}

	$fi = 0.5 * ($fi0 + $fi1); $r = $a * $fi;
	$s = "$fi $r";
	return $s;
}


sub dist {
		my($fi, $a) = @_;
		my($t);

		$t = sqrt($fi*$fi + 1);
		$t = $a * 0.5 * (log($fi + $t) + $fi * $t);

		return $t;
}

sub get_step{
	my($L, $R, $wi) = @_;
	my($a0, $a1, $x0, $x1, $x2, $x3, $y0, $y1, $y2, $y3, $a, $del);

	$R -= $wi; $a0 = 1; $a1 = $R;

	while($a1 - $a0 > 0.001){
		$del = $a1 - $a0; $x0 = $a0; $x1 = $x0 + $del/3; 
		$x2 = $x1 + $del/3; $x3 = $a1;
		$y0 = dist1($x0, $R); $y1 = dist1($x1, $R); 
		$y2 = dist1($x2, $R); $y3 = dist1($x3, $R);
		if($L > $y1){
			$a1 = $x1; 
		}
		elsif($L > $y2){
			$a1 = $x2; $a0 = $x1;
		}
		else{
			$a1 = $x3; $a0 = $x2;
		}
	}

	$a = 0.5 * ($a0 + $a1); 
	return $a;
}

sub dist1{
		my($a, $R) = @_;
		my($t, $r);

		$r = $R/$a; $t = sqrt($r *$r +1);
		$t = 0.5 * $a * (log($r + $t) + $r * $t);

		return $t;
}

sub fromHexToPScolor
{	# GB
	my $hex  = shift;
	$hex =~ s/^#//;
	@hex = (); $ps = "";
	if ($hex =~ /(..)(..)(..)/)
	{
		push @hex, $1, $2, $3;
	}
	
	@nums = (0..9,'A'..'F');
	%nums = map { $nums[$_] => $_ } 0..$#nums;

	foreach $hexThird (@hex)
	{
		$number = 0;
		for( $hexThird =~ /./g )
		{
			$number *= 16;
			$number += $nums{$_};
		}
		$decimal = sprintf("%.2f", $number / 255);
		if ($decimal == 0) { $decimal = 0; }
		if ($decimal == 1) { $decimal = 1; }
		$ps .= "$decimal ";
	}
	$ps =~ s/ $//;	# Drop the last space
	# $clr{$hex} = $ps;
	return $ps;
}