614 lines
12 KiB
C
614 lines
12 KiB
C
#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;
|
|
}
|