#!/usr/bin/perl
#
# Produces list of pixels describing the locations of satellite trails
# and other linear defects in images.
#
# Matthew Hunt <mph@astro.caltech.edu>
#
#	$Id: sat-b-gon.pl,v 1.4 2001/04/29 20:45:18 mph Exp $
#
# Implementation notes:
#	We represent the bad pixel list using a hash called %badpix.  To
#	add a pixel (e.g. 42,69) to the list, we form a scalar key (e.g.
#	"42,69" and set $badpix{"42,69"}.  The value it is set to is
#	irrelevant; we're just taking advantage of Perl's bookkeeping of
#	which keys exist.  We extract the bad pixel list using Perl's
#	"keys" operator which produces a list of all keys.
#
#	The pixel list is grown ("-g" parameter) by iterating over pixels
#	and adding all 8 neighbors of each pixel to the list.

use Getopt::Std;
use POSIX;

getopts('g:hv');			# set $opt_g to -g parameter, etc.

if ($opt_h) {
  print <<EOF;
sat-b-gon [-h] [-g num] endptfile > pixelfile

	-h	This help text
	-g num	Grow trail by radius num pixels
	-v	Verbose

endptfile (one or more, or standard input) consists of one-line descriptions
of each trail to remove.  It should be of the form

x1 y1 x2 y2

Where the trail runs from (x1,y1) to (x2,y2).  A list of bad pixels will
be sent to standard output.  The format is

x y

and is suitable for input to the IRAF "text2mask" procedure.
EOF
  exit 0;
}

$lineno = 0;

# Loop over all specified input
while (<>) {
  ++$lineno;
  next if /^\s*#/;		# skip over comments
  ($x1, $y1, $x2, $y2) = split;
  $ok = 1;
  $ok = 0 if ($x1 == 0 || $y1 == 0 || $x2 == 0 || $y2 == 0);
  warn "bad input line $lineno: $_" unless $ok;
  next unless $ok;

  $m = ($y2-$y1) / ($x2-$x1);	# slope

  if ($opt_v) {
    print STDERR "$x1 $y1 $x2 $y2 m=$m\n";
  }

  if (abs($m) < 1.0) {
    # The "horizontal" case
    # We'll step through every pixel horizontally.
    if ($x2 < $x1) {		# we want increasing x
      $t=$x2; $x2=$x1; $x1=$t;
      $t=$y2; $y2=$y1; $y1=$t;
    }
    $m = ($y2-$y1) / ($x2-$x1);
    $x1 = floor($x1);
    $x2 = ceil($x2);

    $y = $y1;
    for ($x = $x1; $x <= $x2; $x++) {
      $yapprox = floor($y+0.5);	# round y to nearest value
      $idx = sprintf "%d,%d", $x, $yapprox;
      $badpix{$idx} = 1;
      $y += $m;
    }
  } else {
    # The "vertical case"
    if ($y2 < $y1) {
      $t=$x2; $x2=$x1; $x1=$t;
      $t=$y2; $y2=$y1; $y1=$t;
    }
    $m = ($y2-$y1) / ($x2-$x1);
    $y1 = floor($y1);
    $y2 = ceil($y2);

    $x = $x1;
    for ($y = $y1; $y <= $y2; $y++) {
      $xapprox = floor($x+0.5);
      $idx = sprintf "%d,%d", $xapprox, $y;
      $badpix{$idx} = 1;
      $x += 1.0/$m;
    }
  }
}

# Growth

for ($griter = 0; $griter < $opt_g; $griter++) {
  print STDERR "Growth iteration ", $griter+1," /$opt_g\n" if ($opt_v);
  for $groworigin (keys %badpix) {
    ($xor,$yor) = split /,/, $groworigin;
    for ($i = $xor-1; $i <= $xor+1; $i++) {
      for ($j = $yor-1; $j <= $yor+1; $j++) {
	$idx = sprintf "%d,%d", $i, $j;
	$badpix{$idx} = 1;
      }
    }
  }
}

# Emit the output

for $coord (keys %badpix) {
    ($x,$y) = split /,/, $coord;
    printf "%d\t%d\n", $x, $y;
}