#!/usr/local/bin/perl -w
# Ashish Mahabal aam@astro.caltech.edu
# 20 July 2007

# Takes an image (white stars), puts a white circle of specified width at suggested location(s),
# inverts and writes out
# USAGE: $0 -in foo.jpg -out out.jpg [options]

# Uses cimstat and cimstat_mask

BEGIN{ unshift(@INC,"/home/avyakta/aam/usr/local/lib/perl5/site_perl"); }

use Image::Magick;
use strict;
use Carp;
use Getopt::Euclid;
use PDL;

#########################################
# Read in the arguments. No booleans currently.
my $inimage=$ARGV{-inimage};
my $outimage=$ARGV{-outimage};
my $informat=$ARGV{-informat};
my $outformat=$ARGV{-outformat};
my $id=$ARGV{-id};
my $width=$ARGV{-width};
my $ra=$ARGV{-ra};
my $dec=$ARGV{-dec};
my $xcen=$ARGV{-xcen};
my $ycen=$ARGV{-ycen};
my $pointsize=$ARGV{-pointsize};
my $coords=$ARGV{-coords};
my $autoid=$ARGV{-autoid};
my $startid=$ARGV{-startid};
my $separator=$ARGV{-separator};
my $xcol=$ARGV{-xcol};
my $ycol=$ARGV{-ycol};
my $idcol=$ARGV{-idcol};
my $nocen=$ARGV{-nocen};

# Read the input image
my $in1= Image::Magick->new;
if(-e $inimage){ $in1->Read("${informat}:$inimage");}
else{croak "Can not read $inimage\n";}

# Get the image x/y size
my ($xsize,$ysize)=$in1->Get('width','height');

# Is it possible to get the magick without reading the image first?
#$informat=$in1->Get('magick');	# This returns FITS for a fits image

# Make circles thicker if equalizing input image (because its fits)
# May want to switch to the mode Roy uses
# This part is currently still shaky

my $strokewidth=1;
if($informat eq "fits"){
	my $maskimage=$inimage;
	$maskimage=~s/.fits$/.weight.fits/;
	my $stats2=`/home/avyakta/aam/bin/cimstat_mask $inimage $maskimage`;
	my @stats2=split(/\s+/,$stats2);
	my $stats=`/home/donalek/newpipe/bin/cimstat $inimage`;
	my @stats=split(/\s+/,$stats);
	my $median=$stats[15];
	my $stddev=$stats[16];
	my $BlackThreshold=$median-2*$stddev;
	my $WhiteThreshold=$median+3.5*$stddev;
	my @statistics;
	$in1->BlackThreshold(threshold=>$BlackThreshold);
	$in1->WhiteThreshold(threshold=>$WhiteThreshold);
	@statistics=$in1->Statistics();
	$strokewidth=2;

	# 2008-03-01
	system("pdlfits2fits.pl $inimage $maskimage 1");
	my $in2= Image::Magick->new;
	$in2->Read("${informat}:image.fits");
	$in1 = $in2;
}

# Use the central coordinate for id center if none specified
# Use the edges if out of range
if($xcen == 0){$xcen=$xsize/2;}
elsif($xcen < 0){$xcen=-1;}
elsif($xcen > $xsize){$xcen=$xsize+1;}

if($ycen == 0){$ycen=$ysize/2;}
elsif($ycen < 0){$ycen=-1;}
elsif($ycen > $ysize){$ycen=$ysize+1;}

# If ra/dec provided, they get precedence over xcen/ycen
if($informat eq "fits" and $ra != -1 and $dec != -91){
	my $xyarray=`sky2xy $inimage $ra $dec`;
	my @xyarray=split(/\s+/,$xyarray);
	$xcen=$xyarray[4];
	$ycen=$xyarray[5];
}

# Initialize x, y, id arrays to the specified xcen/ycen unless nocen is set

my (@x, @y, @id);
my $idcount=0;
unless($nocen==1){
	$x[$idcount]=$xcen;
	$y[$idcount]=$ycen;
	$id[$idcount]=$id;
	++$idcount;
	}

# If a coords file is provided, attempt to open it
# For lines that are not comments and have at least two columns separated by
# $separator, use them as x and y.
# If xcol, ycol, idcol are specified, ensure that the line has at least that
# many columns and then use those columns instead. NO typechecking done at
# this point. If the columns are non-numeric, that is the users problem

if($coords ne "empty"){
	open my $COORDS,'<',$coords or croak "Can not open $coords";
	while(my $xyid=<$COORDS>){
		$xyid=~s/^\s+//g;		# get rid of leading spaces
		if($xyid=~m/^#/){next;}			# ignore comments
#		print "line  is $xyid\n";
		my @xyid=split(/$separator/,$xyid);	# split in to columns
		if($#xyid<1){next;}		# ignore single column lines
		elsif($#xyid>=$xcol-1 and $#xyid>=$ycol-1 and $#xyid>=$idcol-1){
			#	y coordinates are inverted for fits and other image types
			($x[$idcount],$y[$idcount],$id[$idcount])=($xyid[$xcol],$ysize-$xyid[$ycol]+1,$xyid[$idcol]);
			++$idcount;
			}
		else{ ($x[$idcount],$y[$idcount])=@xyid;$id[$idcount]="";++$idcount; }
		}
	}

my $circle= Image::Magick->new;			# Make a new instance
$circle->Set(size=>"${xsize}x${ysize}");	# Same size as input
$circle->Read('xc:none');			# What is xc exactly?
$circle->Set(matte=>'white');			# Set alpha channel to white

# Once imagemagick is updated from 6.2.4 to 6.3.4 (or beyond) the fill=>'#000000FF'
# in the call to draw circle should be replaced by fill=>'matte'. It is for
# the correct transparency.

#print "idcount is $idcount\n";
for my $i (0..$idcount-1){	# Draw a circle for each x,y pair
	my $cedge=$y[$i]+$width;
	$circle->Draw(primitive=>"circle",points=>"$x[$i],$y[$i] $x[$i],$cedge",
		stroke=>"black",strokewidth=>$strokewidth,fill=>'matte');
#		stroke=>"black",strokewidth=>$strokewidth,fill=>'#000000FF');
	}
$in1->Negate();			# negate the image. Now the stars are black
$in1->Composite(image=>$circle,compose=>'over');	# Add the circles

# If autoid is set (>0) use sequential ids starting from startid (default 0)
if($autoid>0){ for my $i (0..$idcount-1){ $id[$i]=$startid++; } }

for my $i (0..$idcount-1){	# Annotate with ids.
	my $xpos=$x[$i]+$width+1;
	my $ypos=$y[$i]+$width+1;
	$in1->Annotate ( text => $id[$i], pointsize => $pointsize, fill => 'black', x => $xpos, y => $ypos);
	}
$in1->Set(quality=>100);	# Set the quality for writing
$in1->Write("${outformat}:$outimage");	# Write it out
#$circle->Write("${outformat}:circle.jpg");	# Write it out
##############################################
__END__

=head1 NAME

overlay_circle.pl - overlay circle(s) on an image with optional ids

=head1 VERSION

version 0.1

=head1 USAGE

overlay_circle.pl -inimage <INIMAGE> -outimage <OUTIMAGE> [options]

=head1 REQUIRED ARGUMENTS

=over

=item -in[image] [=] <INIMAGE>

Specify input image name 

=for Euclid:
        INIMAGE.type:        string

=item -out[image] [=] <OUTIMAGE>

Specify output image name 

=for Euclid:
        OUTIMAGE.type:        string

=back

=head1 OPTIONS

=over

=item -in[format] [=] <informat>

Specify input image format. Default jpg

=for Euclid:
        informat.type:        string
        informat.default:       'jpg' 

=item -out[format] [=] <outformat>

Specify output image format. Default jpg
Des NOT change the output image extension automagically

=for Euclid:
        outformat.type:        string
        outformat.default:       'jpg' 

=item -coords [=] <COORDS>

Specify file name with coordinates 

=for Euclid:
        COORDS.type:        string
        COORDS.default:       "empty" 

=item -separator [=] <separator>

Specify column separator for coordinate file
Use '\s+' for whitespace

=for Euclid:
        separator.type:        string
        separator.default:     ','

=item -ra [=] <ra> [number]

RA. Used only if inimage is fits.

=for Euclid:
        ra.type:        number >= 0
        ra.type:        number <= 360
	ra.default:    -1 

=item -[dec] [=] <dec> [number]

Dec. Used only if inimage is fits.

=for Euclid:
        dec.type:        number >=-90
        dec.type:        number <=90
	dec.default:    -91 

=item -x[cen] [=] <xcen> [number]

X Position. Defaults to the center.
At the edge if out of range provided.

=for Euclid:
        xcen.type:        number
	xcen.default:    0 

=item -y[cen] [=] <ycen> [number]

Y Position. Defaults to the center.
At the edge if out of range provided.

=for Euclid:
        ycen.type:        number
	ycen.default:    0 

=item -xcol [=] <xcol> [integer]

Xcol in input coords file. Default is 0

=for Euclid:
        xcol.type:        integer >= 0
	xcol.default:     0

=item -ycol [=] <ycol> [integer]

Ycol in input coords file. Default is 1

=for Euclid:
        ycol.type:        integer >= 0
	ycol.default:     1

=item -autoid [=] <autoid> [integer]

number the ids consecutively starting from startid
startid defaults to 0. autoid is 0 (off) by default

=for Euclid:
        autoid.type:        integer >= 0
	autoid.default:     0

=item -startid [=] <startid> [integer]

If autoid is on, the starting pointing of id numbering.
Defaults to 0.

=for Euclid:
        startid.type:        integer
	startid.default:     0

=item -idcol [=] <idcol> [integer]

idcol in input coords file. Default is 2

=for Euclid:
        idcol.type:        integer >= 0
	idcol.default:     2

=item -nocen [=] <nocen> [number]

Whether there should NOT be a circle in the center. Default 0.
i.e. yes, there should be. ID is zero by default

=for Euclid:
        nocen.type:        integer >= 0
	nocen.default:     0

=item -id [=] <id> [string]

ID of the object. Defaults to 0.

=for Euclid:
        id.type:        string
	id.default:     0

=item -p[ointsize] [=] <pointsize> [integer]

Pointsize for font of id. default 20. 0 won't print the id.

=for Euclid:
        pointsize.type:        integer >=0
	pointsize.default:     20

=item -w[idth] [=] <width> [number]

Width of circle

=for Euclid:
        width.type:        number >= 0
	width.default:     10

=item -v

=item --verbose

Print all warnings

=item --version

=item --usage

=item --help

=item --man

Print the usual program information

=back

=head1 AUTHOR

Ashish Mahabal <aam@astro.caltech.edu>

=head1 BUGS

There are undoubtedly serious bugs lurking somewhere in this code.
Bug reports and other feedback are most welcome.

=head1 COPYRIGHT

Copyright (c) 2007, Ashish Mahabal. All Rights Reserved.
This module is free software. It may be used, redistributed
and/or modified under the terms of the Perl Artistic License

(see http://www.perl.com/perl/misc/Artistic.html)