spacecruft
/
sattools
1
0
Fork 0
Code Releases Activity
sattools/src/pgm2fits.c

682 lines
14 KiB
C
Raw Permalink Blame History

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "qfits.h"
#include <getopt.h>
struct image
{
int nx, ny;
char timestamp[64];
unsigned char *c;
};
struct fourframe
{
int nx, ny, nt, nlayer;
char timestamp[64], observer[64];
double mjd;
float *z, *dt;
int cospar, type;
};
int fgetline (FILE * file, char *s, int lim);
void write_fits (char *filename, struct fourframe ff);
struct image read_pgm (char *filename, int *status);
double nfd2mjd (char *date);
double date2mjd (int year, int month, double day);
void write_pgm (char *filename, struct fourframe ff);
void mjd2date (double mjd, char *date);
void
usage (void)
{
printf ("pgm2fits p:w:h:s:n:Dd:x:y:c:o:gm:t:r:I\n\n");
printf ("-p image prefix\n");
printf ("-w image width in pixels\n");
printf ("-h image height in pixels\n");
printf ("-s number of first image to process\n");
printf ("-n number of images to process\n");
printf ("-D toggle for creating dark frame\n");
printf ("-d toggle for creating dark frame\n");
printf ("-d filename of dark frame to substract\n");
printf ("-m filename of mask frame to apply\n");
printf ("-x tracking rate in x (pix/s)\n");
printf ("-y tracking rate in y (pix/s)\n");
printf ("-c COSPAR [default 4553]\n");
printf ("-o observer [default \"Cees Bassa\"]\n");
printf ("-g toggle for guiding??\n");
printf ("-t time stamp of first image [YYYY-MM-DDTHH:MM:SS.SSS]\n");
printf ("-r frame rate (frames/s)\n");
printf ("-I integer output\n");
exit (0);
}
int
main (int argc, char *argv[])
{
int i, j, k, l, m, k0 = 1, status, nt, darkout = 0, darkin = 0, track =
0, maskin = 0;
int n, n0, di, dj, npix;
struct image *img, drk, msk;
struct fourframe ff;
char filename[128], nfd[64];
float s1, s2, z;
float avg, std, max, cnt, *trk;
int *wt;
int arg = 0;
char *path, *darkfile, *maskfile;
double mjd, mjd0 = 0.0;
float dxdn = 0.0, dydn = 0.0, dx, dy;
int guide = 0, timereset = 0;
float framerate = 25.0;
char *env;
// Set defaults
env = getenv ("ST_COSPAR");
ff.cospar = atoi (env);
strcpy (ff.observer, "Cees Bassa");
ff.type = 0;
// Decode options
if (argc > 1)
{
while ((arg =
getopt (argc, argv, "p:w:h:s:n:Dd:x:y:c:o:gm:t:r:I")) != -1)
{
switch (arg)
{
case 'p':
path = optarg;
break;
case 'g':
guide = 1;
break;
case 'w':
ff.nx = atoi (optarg);
break;
case 'h':
ff.ny = atoi (optarg);
break;
case 'I':
ff.type = 1;
break;
case 'c':
ff.cospar = atoi (optarg);
break;
case 'o':
strcpy (ff.observer, optarg);
break;
case 's':
k0 = atoi (optarg);
break;
case 'n':
nt = atoi (optarg);
break;
case 'D':
darkout = 1;
break;
case 'd':
darkin = 1;
darkfile = optarg;
break;
case 'm':
maskin = 1;
maskfile = optarg;
break;
case 'x':
dxdn = atof (optarg);
track = 1;
break;
case 'y':
dydn = atof (optarg);
track = 1;
break;
case 't':
strcpy (nfd, optarg);
mjd0 = nfd2mjd (nfd);
timereset = 1;
break;
case 'r':
framerate = atof (optarg);
timereset = 1;
break;
default:
usage ();
}
}
}
else
{
usage ();
}
// Add layer
if (track == 1)
ff.nlayer = 5;
else
ff.nlayer = 4;
// Allocate
ff.z = (float *) malloc (ff.nlayer * sizeof (float) * ff.nx * ff.ny);
ff.dt = (float *) malloc (sizeof (float) * nt);
trk = (float *) malloc (sizeof (float) * ff.nx * ff.ny);
wt = (int *) malloc (sizeof (float) * ff.nx * ff.ny);
img = (struct image *) malloc (sizeof (struct image) * nt);
// Read dark file
if (darkin == 1)
drk = read_pgm (darkfile, &status);
// Read mask file
if (maskin == 1)
msk = read_pgm (maskfile, &status);
// Loop over files
for (k = 0, l = 0; k < nt; k++)
{
sprintf (filename, "%s%06d.pgm", path, k + k0);
img[l] = read_pgm (filename, &status);
// Reset time
if (timereset == 1)
{
mjd = mjd0 + (double) (k + k0) / (86400.0 * framerate);
mjd2date (mjd, img[l].timestamp);
}
ff.dt[l] =
86400.0 * (nfd2mjd (img[l].timestamp) - nfd2mjd (img[0].timestamp));
if (status == 0)
{
printf ("Read %s\n", filename);
l++;
}
else
{
break;
}
}
ff.nt = l;
strcpy (ff.timestamp, img[0].timestamp);
ff.mjd = nfd2mjd (img[0].timestamp);
printf ("Accumulating image statistics\n");
// Loop over pixels
for (i = 0; i < ff.nx; i++)
{
for (j = 0; j < ff.ny; j++)
{
n = i + ff.nx * j;
s1 = 0.0;
s2 = 0.0;
max = 0.0;
cnt = 0.0;
// Loop over images
for (k = 0; k < ff.nt; k++)
{
if (darkin == 0)
z = (float) img[k].c[n];
else if (darkin == 1)
z = (float) (img[k].c[n] - drk.c[n]);
s1 += z;
s2 += z * z;
if (z > max)
{
max = z;
cnt = (float) k;
}
}
s1 -= max;
s2 -= max * max;
avg = s1 / (float) (ff.nt - 1);
std = sqrt ((s2 - s1 * avg) / (float) (ff.nt - 2));
// Reset masked pixels
if (maskin == 1 && msk.c[n] == 0.0)
{
avg = 0.0;
std = 0.0;
max = 0.0;
cnt = 128.0;
}
for (m = 0; m < ff.nlayer; m++)
{
l = i + (ff.ny - j - 1) * ff.nx + m * ff.nx * ff.ny;
if (m == 0)
ff.z[l] = avg;
if (m == 1)
ff.z[l] = std;
if (m == 2)
ff.z[l] = max;
if (m == 3)
ff.z[l] = cnt;
if (m == 4)
ff.z[l] = avg;
}
}
}
// Create tracked layer
if (track == 1)
{
printf ("Creating tracked layer\n");
// Set weights
for (i = 0; i < ff.nx * ff.ny; i++)
{
wt[i] = 0;
trk[i] = 0.0;
}
// Loop over frames
for (l = 0; l < ff.nt; l++)
{
// Offset
dx = dxdn * (l - ff.nt / 2);
dy = dydn * (l - ff.nt / 2);
// Integer offset
di = (int) floor (dx + 0.5);
dj = (int) floor (dy + 0.5);
// Set
for (i = 0; i < ff.nx; i++)
{
for (j = 0; j < ff.ny; j++)
{
k = i + ff.nx * j;
k0 = i + di + ff.nx * (j + dj);
if (i + di > 0 && i + di < ff.nx && j + dj > 0
&& j + dj < ff.ny)
{
wt[k] += 1;
trk[k] += (float) img[l].c[k0];
}
}
}
}
// Save layer
for (i = 0; i < ff.nx; i++)
{
for (j = 0; j < ff.ny; j++)
{
k = i + ff.nx * j;
if (guide == 0)
l = i + (ff.ny - j - 1) * ff.nx;
else
l = i + (ff.ny - j - 1) * ff.nx + 4 * ff.nx * ff.ny;
if (wt[k] > 0)
ff.z[l] = trk[k] / (float) wt[k];
else
ff.z[l] = trk[k];
}
}
}
// Write fits
if (ff.timestamp != NULL)
sprintf (filename, "%s.fits", ff.timestamp);
else
strcpy (filename, "test.fits");
write_fits (filename, ff);
// Write dark frame
if (darkout == 1)
write_pgm ("dark.pgm", ff);
// Free
for (k = 0; k < ff.nt; k++)
free (img[k].c);
free (trk);
free (img);
free (wt);
free (ff.dt);
free (ff.z);
return 0;
}
// Read pgm file
struct image
read_pgm (char *filename, int *status)
{
int i;
struct image img;
FILE *file;
char hbuf[64];
// Open file
file = fopen (filename, "rb");
if (file == NULL)
{
*status = 1;
return img;
}
// Read PGM format
fgetline (file, hbuf, 64);
if (strcmp (hbuf, "P5") != 0)
{
printf ("Not a valid PGM file!\n");
exit (0);
}
// Read timestamp/image size
fgetline (file, hbuf, 64);
if (strstr (hbuf, "#") != NULL)
{
strcpy (img.timestamp, hbuf + 2);
fgetline (file, hbuf, 64);
sscanf (hbuf, "%d %d", &img.nx, &img.ny);
}
else
{
strcpy (img.timestamp, "2012-01-01T00:00:00");
sscanf (hbuf, "%d %d", &img.nx, &img.ny);
}
fgetline (file, hbuf, 64);
// Allocate
img.c = (unsigned char *) malloc (sizeof (unsigned char) * img.nx * img.ny);
// Read
fread (img.c, 1, img.nx * img.ny, file);
// Close file
fclose (file);
*status = 0;
return img;
}
// Get line
int
fgetline (FILE * file, char *s, int lim)
{
int c, i = 0;
while (--lim > 0 && (c = fgetc (file)) != EOF && c != '\n')
s[i++] = c;
// if (c=='\n')
// s[i++]=c;
s[i] = '\0';
return i;
}
// Write fits file
void
write_fits (char *filename, struct fourframe ff)
{
int i, j, k, l;
float *fbuf;
int *ibuf;
qfitsdumper qd;
qfits_header *qh;
char key[FITS_LINESZ + 1];
char val[FITS_LINESZ + 1];
char com[FITS_LINESZ + 1];
char lin[FITS_LINESZ + 1];
FILE *file;
// Create FITS header
qh = qfits_header_default ();
// Add stuff
if (ff.type == 1)
qfits_header_add (qh, "BITPIX", "8", " ", NULL);
else
qfits_header_add (qh, "BITPIX", "-32", " ", NULL);
qfits_header_add (qh, "NAXIS", "3", " ", NULL);
sprintf (val, "%i", ff.nx);
qfits_header_add (qh, "NAXIS1", val, " ", NULL);
sprintf (val, "%i", ff.ny);
qfits_header_add (qh, "NAXIS2", val, " ", NULL);
sprintf (val, "%i", ff.nlayer);
qfits_header_add (qh, "NAXIS3", val, " ", NULL);
qfits_header_add (qh, "BSCALE", "1.0", " ", NULL);
qfits_header_add (qh, "BZERO", "0.0", " ", NULL);
qfits_header_add (qh, "DATAMAX", "255.0", " ", NULL);
qfits_header_add (qh, "DATAMIN", "0.0", " ", NULL);
sprintf (val, "%s", ff.timestamp);
qfits_header_add (qh, "DATE-OBS", val, " ", NULL);
// MJD-OBS
sprintf (val, "%lf", ff.mjd);
qfits_header_add (qh, "MJD-OBS", val, " ", NULL);
sprintf (val, "%f", ff.dt[ff.nt - 1] - ff.dt[0]);
qfits_header_add (qh, "EXPTIME", val, " ", NULL);
sprintf (val, "%d", ff.nt);
qfits_header_add (qh, "NFRAMES", val, " ", NULL);
// Astrometry keywors
sprintf (val, "%f", ff.nx / 2.0);
qfits_header_add (qh, "CRPIX1", val, " ", NULL);
sprintf (val, "%f", ff.ny / 2.0);
qfits_header_add (qh, "CRPIX2", val, " ", NULL);
qfits_header_add (qh, "CRVAL1", "0.0", " ", NULL);
qfits_header_add (qh, "CRVAL2", "0.0", " ", NULL);
qfits_header_add (qh, "CD1_1", "0.0", " ", NULL);
qfits_header_add (qh, "CD1_2", "0.0", " ", NULL);
qfits_header_add (qh, "CD2_1", "0.0", " ", NULL);
qfits_header_add (qh, "CD2_2", "0.0", " ", NULL);
qfits_header_add (qh, "CTYPE1", "'RA---TAN'", " ", NULL);
qfits_header_add (qh, "CTYPE2", "'DEC--TAN'", " ", NULL);
qfits_header_add (qh, "CUNIT1", "'deg'", " ", NULL);
qfits_header_add (qh, "CUNIT2", "'deg'", " ", NULL);
qfits_header_add (qh, "CRRES1", "0.0", " ", NULL);
qfits_header_add (qh, "CRRES2", "0.0", " ", NULL);
qfits_header_add (qh, "EQUINOX", "2000.0", " ", NULL);
qfits_header_add (qh, "RADECSYS", "ICRS", " ", NULL);
sprintf (val, "%d", ff.cospar);
qfits_header_add (qh, "COSPAR", val, " ", NULL);
sprintf (val, "'%s'", ff.observer);
qfits_header_add (qh, "OBSERVER", val, " ", NULL);
// Add timestamps
for (k = 0; k < ff.nt; k++)
{
sprintf (key, "DT%04d", k);
sprintf (val, "%f", ff.dt[k]);
qfits_header_add (qh, key, val, " ", NULL);
}
// Dummy keywords
for (k = 0; k < 10; k++)
{
sprintf (key, "DUMY%03d", k);
qfits_header_add (qh, key, "0.0", " ", NULL);
}
// Dump fitsheader
// qfits_header_dump(qh,stdout);
// Dump to file
file = fopen (filename, "w");
qfits_header_dump (qh, file);
fclose (file);
// Fill buffer
fbuf = malloc (ff.nlayer * ff.nx * ff.ny * sizeof (float));
ibuf = malloc (ff.nlayer * ff.nx * ff.ny * sizeof (int));
for (i = 0, l = 0; i < ff.nx; i++)
{
for (j = ff.ny - 1; j >= 0; j--)
{
for (k = 0; k < ff.nlayer; k++)
{
fbuf[l] = ff.z[l];
ibuf[l] = (int) ff.z[l];
l++;
}
}
}
// Set parameters
qd.filename = filename;
qd.npix = ff.nlayer * ff.nx * ff.ny;
if (ff.type == 1)
{
qd.ptype = PTYPE_INT;
qd.ibuf = ibuf;
qd.out_ptype = BPP_8_UNSIGNED;
}
else
{
qd.ptype = PTYPE_FLOAT;
qd.fbuf = fbuf;
qd.out_ptype = -32;
}
// Dump
qfits_pixdump (&qd);
free (fbuf);
free (ibuf);
return;
}
// Compute Julian Day from Date
double
date2mjd (int year, int month, double day)
{
int a, b;
double jd;
if (month < 3)
{
year--;
month += 12;
}
a = floor (year / 100.);
b = 2. - a + floor (a / 4.);
if (year < 1582)
b = 0;
if (year == 1582 && month < 10)
b = 0;
if (year == 1582 && month == 10 && day <= 4)
b = 0;
jd =
floor (365.25 * (year + 4716)) + floor (30.6001 * (month + 1)) + day + b -
1524.5;
return jd - 2400000.5;
}
// nfd2mjd
double
nfd2mjd (char *date)
{
int year, month, day, hour, min;
double mjd, dday;
float sec;
sscanf (date, "%04d-%02d-%02dT%02d:%02d:%f", &year, &month, &day, &hour,
&min, &sec);
dday = day + hour / 24.0 + min / 1440.0 + sec / 86400.0;
mjd = date2mjd (year, month, dday);
return mjd;
}
// Write pgm file
void
write_pgm (char *filename, struct fourframe ff)
{
int i, j, k;
FILE *file;
file = fopen (filename, "w");
fprintf (file, "P5\n# 2013-01-01T00:00:00\n%d %d\n255\n", ff.nx, ff.ny);
for (j = 0; j < ff.ny; j++)
{
for (i = 0; i < ff.nx; i++)
{
k = i + (ff.ny - j - 1) * ff.nx;
fprintf (file, "%c", (char) ff.z[k]);
}
}
fclose (file);
return;
}
// Compute Date from Julian Day
void
mjd2date (double mjd, char *date)
{
double f, jd, dday;
int z, alpha, a, b, c, d, e;
int year, month, day, hour, min;
double sec, x;
jd = mjd + 2400000.5;
jd += 0.5;
z = floor (jd);
f = fmod (jd, 1.);
if (z < 2299161)
a = z;
else
{
alpha = floor ((z - 1867216.25) / 36524.25);
a = z + 1 + alpha - floor (alpha / 4.);
}
b = a + 1524;
c = floor ((b - 122.1) / 365.25);
d = floor (365.25 * c);
e = floor ((b - d) / 30.6001);
dday = b - d - floor (30.6001 * e) + f;
if (e < 14)
month = e - 1;
else
month = e - 13;
if (month > 2)
year = c - 4716;
else
year = c - 4715;
day = (int) floor (dday);
x = 24.0 * (dday - day);
x = 3600. * fabs (x);
sec = fmod (x, 60.);
x = (x - sec) / 60.;
min = fmod (x, 60.);
x = (x - min) / 60.;
hour = x;
sprintf (date, "%04d-%02d-%02dT%02d:%02d:%06.3f", year, month, day, hour,
min, sec);
return;
}
View Git Blame Copy Permalink