#!/usr/bin/env /proj/axaf/bin/perl
# --8<--8<--8<--8<--
#
# Copyright (C) 2008 Smithsonian Astrophysical Observatory
#
# This file is part of rawmedadd
#
# rawmedadd is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
#
# -->8-->8-->8-->8--
use strict;
use warnings;
use PDL;
use Getopt::Long;
use File::Spec::Functions qw( catfile );
my $VERSION = '1.0.0';
my %opt;
my $have_Math_Prime = eval "use Math::Prime::XS; 1;";
eval { main() };
if ( $@ )
{
warn $@;
exit 1;
}
exit 0;
sub main {
parse_opts();
help(1) if $opt{help};
help(2) if $opt{usage};
do { print "$0 $VERSION\n"; return }
if $opt{version};
# list of images
my @images = sort glob( catfile( $opt{dir}, '*_img0.fits*') );
die( "no images in directory `$opt{dir}'\n" )
unless @images;
my %winfo = get_window_info( $opt{dir} );
if ( $opt{autonmed} )
{
$opt{nmed} ||= 3;
my $nmed = $opt{nmed} ? prime( scalar @images ) : 0;
if ( defined $nmed && $nmed != $opt{nmed} )
{
$opt{nmed} = $nmed;
warn "adjusted nmed to fully populate groups (nmed is now $opt{nmed})\n";
}
elsif ( ! defined $nmed )
{
my $ndiscard = scalar @images % $opt{nmed};
warn "unable to tune nmed to fully populate groups; last incomplete group of $ndiscard will be discarded\n";
}
}
elsif ( $opt{nmed} )
{
if ( $opt{nmed} == -1 )
{
$opt{nmed} = scalar @images if $opt{nmed} == -1;
if ( $opt{nmed} %2 )
{
warn "nmed ($opt{nmed}) is not an odd number; nmed reduced accordingly\n";
$opt{nmed}--;
}
}
elsif ( $opt{nmed} > @images )
{
warn "nmed ($opt{nmed}) > number of images (", scalar @images, "); adjusted accordingly\n";
$opt{nmed} = @images;
}
elsif ( my $ndiscard = @images % $opt{nmed} )
{
warn "nmed of $opt{nmed} will not fully populate groups; last incomplete group of $ndiscard will be discarded\n";
}
}
# output image
my $sum = zeroes( float(), $winfo{nx}, $winfo{ny} );
my $nslices = $opt{nmed} ? int( @images / $opt{nmed} ) : @images;
my $cube = exists $opt{cube} ? zeroes( float(), $winfo{nx}, $winfo{ny}, $nslices ) : undef;
if ( $opt{nmed} )
{
addmedian( $sum, $cube, \%winfo, @images );
}
else
{
addall( $sum, $cube, \%winfo, @images );
}
$sum->wfits( $opt{outfile} );
$cube->wfits( $opt{cube} )
if defined $cube;
}
sub addmedian {
my ( $sum, $cube, $winfo ) = splice( @_, 0, 3);
my $exposure = 0;
my $group = zeroes( float(), $opt{nmed}, $winfo->{nx}, $winfo->{ny} );
my $idx = 0;
while ( @_ )
{
print scalar @_, "\n";
my @images = splice( @_, 0, $opt{nmed} );
last if @images != $opt{nmed};
my @exposures;
# grab images
for my $gidx ( 0..$opt{nmed}-1 )
{
my $img = getimg( $images[$gidx], $winfo );
push @exposures, $img->gethdr->{EXPTIME};
$group->mslice([$gidx],[],[]) .= $img->dummy(0);
}
$exposure += pdl(@exposures)->avg;
# calculate the median of each pixel
my $median = $group->medover;
$sum += $median;
$cube->mslice([],[],[$idx]) .= $median
if defined $cube;
}
continue
{
$idx++;
}
my $hdr = $sum->fhdr;
$hdr->{EXPTIME} = $exposure;
$hdr->{NMED} = $opt{nmed};
$hdr->{NFRAMES} = $idx;
}
sub addall{
my ( $sum, $cube, $winfo ) = splice( @_, 0, 3);
my $exposure = 0;
my $idx = 0;
for my $file ( @_ )
{
print scalar @_ - $idx, "\n";
my $img = getimg( $file, $winfo );
$sum += $img;
$exposure += $img->gethdr->{EXPTIME};
$cube->mslice([],[],[$idx]) .= $img
if defined $cube;
}
continue
{
$idx++;
}
my $hdr = $sum->fhdr;
$hdr->{EXPTIME} = $exposure;
$hdr->{NMED} = $opt{nmed};
$hdr->{NFRAMES} = $idx;
}
sub getimg
{
my ( $file, $winfo ) = @_;
my $fits = rfits( $file, { hdrcpy => 1} )->float;
# bias columns
my $biascols = $fits->mslice( [1024, 1024+63], [] );
for my $node ( 0..3 )
{
my $img = $fits->mslice( [ $node * 256,
$node * 256 + 255], [] );
my $bias = $biascols->mslice( [ $node * 16,
$node * 16 + 15], [] )->float->average;
$img -= $bias->dummy(0);
}
# get windowed image
my $img = $fits->mslice( [ $winfo->{xmin}, $winfo->{xmax} ],
[ $winfo->{ymin}, $winfo->{ymax} ] )->copy ;
my $hdr = $img->gethdr;
$hdr->{EXPTIME} = $hdr->{TSTOP} - $hdr->{TSTART};
delete @{$hdr}{ qw( TSTART TSTOP ) };
return $img;
}
sub get_window_info
{
my ( $dir ) = @_;
# window file
my @win = sort glob( catfile( $dir, '*_win0.fits*') );
die( "no or more than one window file: @win\n" )
if @win != 1;
my $win = shift @win;
my $tbl = rfits( $win . '[1]' );
my %info;
# grab the first row; do some magic to get the indices correct
$info{nx} = $tbl->{CCDCOL}->at(0) + 1;
$info{ny} = $tbl->{CCDROW}->at(0) + 1;
$info{xmin} = $tbl->{LL_CCDX}->at(0);
$info{ymin} = $tbl->{LL_CCDY}->at(0);
$info{xmax} = $info{xmin} + $info{nx} - 1;
$info{ymax} = $info{ymin} + $info{ny} - 1;
return %info;
}
sub parse_opts {
%opt = (
nmed => 5,
outfile => 'sum.fits',
autonmed => 0,
);
eval {
local $SIG{__WARN__} = sub { die $_[0] };
GetOptions( \%opt,
qw/
autonmed
cube=s
nmed=i
outfile=s
help
usage
version
/);
};
die $@ if $@;
return if $opt{help} || $opt{usage} || $opt{version};
die( "number of images in median group must be zero or an odd number greater than one\n" )
if $opt{nmed} != 0 && ( $opt{nmed} == 1 || $opt{nmed}%2 == 0 );
$opt{dir} = shift @ARGV || '.';
die( "extra stuff on the command line: @ARGV\n" )
if @ARGV;
$opt{outfile} .= '.fits' unless $opt{outfile} =~ /[.]fits$/;
}
sub help
{
my ( $verbose ) = @_;
# verbosity = 2 causes pod2usage to call perldoc, which already pages
require Pod::Usage;
Pod::Usage::pod2usage ( { -exitval => 0, -verbose => $verbose } );
}
sub prime
{
my ( $n ) = @_;
my $prime;
if ( $have_Math_Prime )
{
my @primes = Math::Prime::XS::primes( $opt{nmed}, $n );
$prime = shift @primes;
}
else
{
$prime = $opt{nmed};
while( $prime < $n )
{
last if ($n%$prime) == 0;
$prime+=2;
}
$prime = undef if $prime >= $n;
}
return $prime;
}
__END__
=head1 NAME
rawmedadd - sum ACIS RAW frames with median cleanup
=head1 SYNOPSIS
B I F
=head1 DESCRIPTION
B adds ACIS data in C data mode. It can perform
median cleaning, adding the median images.
It operates on an entire directory of images considering files with
the suffix of F<_img0.fits*> to be images. It uses the F<_win0.fits*>
file to determine the size of the data window.
If the B<--nmed> option is non-zero, B will group the images
into groups of size C, create a single median image from each group,
then add the median images.
=head1 ARGUMENTS AND OPTIONS
B uses long options. Options may be abbreviated, and the "="
character shown below in the option templates is optional.
The directory in which the ACIS data is to be found should be the last
argument on the command line. If not specified, the current directory
is used.
The following options are recognized:
=over
=item --autonmed
If true, the number of images in a median group will be automatically
determined, using the value of B<--nmed> as a starting point. It defaults to 0.
=item --cube=F
If specified, a three dimensional FITS image cube will be written to the
named file, containing the stack of median images (or of all images if the
B argument is C<0>.
=item --nmed=I
The number of images in a median group. Defaults to C<5>. If C<0>,
the images are added directly. If C<-1>, the median of all of the
images is used.
=item --outfile=F
The output FITS image name. A suffix of C<.fits> will be automatically supplied.
This defaults to F.
=back
=head1 COPYRIGHT & LICENSE
Copyright 2008 Smithsonian Astrophysical Observatory
This software is released under the GNU General Public License. You
may find a copy at
http://www.fsf.org/copyleft/gpl.html
=head1 AUTHOR
Diab Jerius Edjerius@cfa.harvard.eduE
=cut