326 lines
7.3 KiB
C
326 lines
7.3 KiB
C
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include "qfits.h"
|
|
#include <cpgplot.h>
|
|
#include <wcslib/cel.h>
|
|
#include <jpeglib.h>
|
|
#include <getopt.h>
|
|
|
|
struct image {
|
|
int naxis1,naxis2,naxis3;
|
|
float *z;
|
|
float zmin,zmax;
|
|
double ra0,de0;
|
|
float avg,std;
|
|
float x0,y0;
|
|
float a[3],b[3],xrms,yrms;
|
|
float exptime;
|
|
double mjd;
|
|
char nfd[32];
|
|
int cospar;
|
|
};
|
|
struct jpeg_image {
|
|
int nx,ny,nz;
|
|
float *z;
|
|
};
|
|
struct jpeg_image read_jpg(char *filename);
|
|
void write_jpg(char *filename,struct jpeg_image img);
|
|
struct image read_fits(char *filename,int pnum);
|
|
void forward(double ra0,double de0,double ra,double de,double *x,double *y);
|
|
void reverse(double ra0,double de0,double x,double y,double *ra,double *de);
|
|
|
|
int main(int argc,char *argv[])
|
|
{
|
|
int i,j,k,l,m;
|
|
struct image img;
|
|
struct jpeg_image jpg,out;
|
|
double rx,ry,ra,de,rx0,ry0;
|
|
double x,y,d;
|
|
double drx=-10.0,dry=10.0;
|
|
double ra0=237.0,de0=12.5;
|
|
int arg=0;
|
|
char *fitsfile,*jpgfile,*outfile;
|
|
|
|
// Decode options
|
|
while ((arg=getopt(argc,argv,"j:f:o:R:D:s:"))!=-1) {
|
|
switch(arg) {
|
|
|
|
case 'j':
|
|
jpgfile=optarg;
|
|
break;
|
|
|
|
case 'f':
|
|
fitsfile=optarg;
|
|
break;
|
|
|
|
case 'o':
|
|
outfile=optarg;
|
|
break;
|
|
|
|
case 'R':
|
|
ra0=atof(optarg);
|
|
break;
|
|
|
|
case 'D':
|
|
de0=atof(optarg);
|
|
break;
|
|
|
|
case 's':
|
|
dry=atof(optarg);
|
|
drx=-dry;
|
|
break;
|
|
|
|
default:
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
|
|
// Read image
|
|
img=read_fits(fitsfile,0);
|
|
jpg=read_jpg(jpgfile);
|
|
|
|
out.nx=3000;
|
|
out.ny=6000;
|
|
out.nz=3;
|
|
|
|
/*
|
|
img.x0*=4.0;
|
|
img.y0*=4.0;
|
|
img.a[1]/=4.0;
|
|
img.a[2]/=4.0;
|
|
img.b[1]/=4.0;
|
|
img.b[2]/=4.0;
|
|
*/
|
|
out.z=(float *) malloc(sizeof(float)*out.nx*out.ny*out.nz);
|
|
|
|
for (i=0;i<out.nx;i++) {
|
|
for (j=0;j<out.ny;j++) {
|
|
// Set rx,ry
|
|
rx=drx*(float) (i-0.5*out.nx);
|
|
ry=dry*(float) (j-0.5*out.ny);
|
|
|
|
// Obtain ra/dec for output image
|
|
reverse(ra0,de0,rx,ry,&ra,&de);
|
|
|
|
// Obtain rx/ry for input image
|
|
forward(img.ra0,img.de0,ra,de,&rx0,&ry0);
|
|
|
|
// Compute pixel position
|
|
d=img.a[1]*img.b[2]-img.a[2]*img.b[1];
|
|
x=(+rx0*img.b[2]-ry0*img.a[2])/d+img.x0;
|
|
y=(-rx0*img.b[1]+ry0*img.a[1])/d+img.y0;
|
|
|
|
// Fill image
|
|
for (k=0;k<jpg.nz;k++) {
|
|
l=out.nz*(i+out.nx*(out.ny-j-1))+k;
|
|
m=jpg.nz*((int) x+jpg.nx*(int) (jpg.ny-y-1))+k;
|
|
if (x>0.0 && x<jpg.nx && y>0.0 && y<jpg.ny)
|
|
out.z[l]=jpg.z[m];
|
|
else
|
|
out.z[l]=0.0;
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
// Dump
|
|
write_jpg(outfile,out);
|
|
|
|
// Free
|
|
free(img.z);
|
|
free(jpg.z);
|
|
free(out.z);
|
|
|
|
return 0;
|
|
}
|
|
|
|
// Read fits image
|
|
struct image read_fits(char *filename,int pnum)
|
|
{
|
|
int i,j,k,l,m;
|
|
qfitsloader ql;
|
|
char key[FITS_LINESZ+1] ;
|
|
struct image img;
|
|
float s1,s2,avg,std;
|
|
|
|
// Set plane
|
|
ql.xtnum = 0;
|
|
ql.pnum = pnum;
|
|
|
|
// Set loadtype
|
|
ql.ptype = PTYPE_FLOAT;
|
|
|
|
// Set filename
|
|
ql.filename=filename;
|
|
|
|
// Image size
|
|
img.naxis1=atoi(qfits_query_hdr(filename,"NAXIS1"));
|
|
img.naxis2=atoi(qfits_query_hdr(filename,"NAXIS2"));
|
|
|
|
// MJD
|
|
img.mjd=atof(qfits_query_hdr(filename,"MJD-OBS"));
|
|
strcpy(img.nfd,qfits_query_hdr(filename,"DATE-OBS"));
|
|
img.exptime=atof(qfits_query_hdr(filename,"EXPTIME"));
|
|
|
|
// COSPAR ID
|
|
img.cospar=atoi(qfits_query_hdr(filename,"COSPAR"));
|
|
|
|
// Transformation
|
|
img.mjd=atof(qfits_query_hdr(filename,"MJD-OBS"));
|
|
img.ra0=atof(qfits_query_hdr(filename,"CRVAL1"));
|
|
img.de0=atof(qfits_query_hdr(filename,"CRVAL2"));
|
|
img.x0=atof(qfits_query_hdr(filename,"CRPIX1"));
|
|
img.y0=atof(qfits_query_hdr(filename,"CRPIX2"));
|
|
img.a[0]=0.0;
|
|
img.a[1]=3600.0*atof(qfits_query_hdr(filename,"CD1_1"));
|
|
img.a[2]=3600.0*atof(qfits_query_hdr(filename,"CD1_2"));
|
|
img.b[0]=0.0;
|
|
img.b[1]=3600.0*atof(qfits_query_hdr(filename,"CD2_1"));
|
|
img.b[2]=3600.0*atof(qfits_query_hdr(filename,"CD2_2"));
|
|
img.xrms=3600.0*atof(qfits_query_hdr(filename,"CRRES1"));
|
|
img.yrms=3600.0*atof(qfits_query_hdr(filename,"CRRES2"));
|
|
|
|
// Initialize load
|
|
if (qfitsloader_init(&ql) != 0)
|
|
printf("Error initializing data loading\n");
|
|
|
|
// Test load
|
|
if (qfits_loadpix(&ql) != 0)
|
|
printf("Error loading actual data\n");
|
|
|
|
// Allocate image memory
|
|
img.z=(float *) malloc(sizeof(float) * img.naxis1*img.naxis2);
|
|
|
|
// Fill z array
|
|
for (i=0,l=0,m=0;i<img.naxis1;i++) {
|
|
for (j=0;j<img.naxis2;j++) {
|
|
img.z[l]=ql.fbuf[l];
|
|
l++;
|
|
}
|
|
}
|
|
|
|
// Get levels
|
|
for (i=0,s1=0.0,s2=0.0;i<img.naxis1*img.naxis2;i++) {
|
|
s1+=img.z[i];
|
|
s2+=img.z[i]*img.z[i];
|
|
}
|
|
img.avg=s1/(float) (img.naxis1*img.naxis2);
|
|
img.std=sqrt(s2/(float) (img.naxis1*img.naxis2)-img.avg*img.avg);
|
|
img.zmin=img.avg-4.0*img.std;
|
|
img.zmax=img.avg+12.0*img.std;
|
|
|
|
return img;
|
|
}
|
|
|
|
struct jpeg_image read_jpg(char *filename)
|
|
{
|
|
int i=0,j,k,l,m;
|
|
unsigned long location=0;
|
|
struct jpeg_image img;
|
|
struct jpeg_decompress_struct cinfo;
|
|
struct jpeg_error_mgr jerr;
|
|
JSAMPROW row_pointer[1];
|
|
unsigned char *raw_image=NULL;
|
|
FILE *file;
|
|
|
|
// Open file
|
|
file=fopen(filename,"rb");
|
|
if (!file)
|
|
perror("Error opening file");
|
|
|
|
// Get header info, decompress
|
|
cinfo.err=jpeg_std_error(&jerr);
|
|
jpeg_create_decompress(&cinfo);
|
|
jpeg_stdio_src(&cinfo,file);
|
|
jpeg_read_header(&cinfo,TRUE);
|
|
jpeg_start_decompress(&cinfo);
|
|
|
|
// Allocate memory
|
|
raw_image=(unsigned char *) malloc(cinfo.output_width*cinfo.output_height*cinfo.num_components);
|
|
|
|
// Read image, one scan at a time
|
|
row_pointer[0]=(unsigned char *) malloc(cinfo.output_width*cinfo.num_components);
|
|
while(cinfo.output_scanline<cinfo.image_height) {
|
|
jpeg_read_scanlines(&cinfo,row_pointer,1);
|
|
for(i=0;i<cinfo.image_width*cinfo.num_components;i++)
|
|
raw_image[location++]=row_pointer[0][i];
|
|
}
|
|
// wrap up decompression, destroy objects, free pointers and close open files
|
|
jpeg_finish_decompress(&cinfo);
|
|
jpeg_destroy_decompress(&cinfo);
|
|
|
|
// Copy image to image struct
|
|
img.nx=cinfo.image_width;
|
|
img.ny=cinfo.image_height;
|
|
img.nz=cinfo.num_components;
|
|
img.z=(float *) malloc(sizeof(float)*img.nx*img.ny*img.nz);
|
|
|
|
// Fill image
|
|
for (i=0;i<img.nx;i++) {
|
|
for (j=0;j<img.ny;j++) {
|
|
for (k=0;k<img.nz;k++) {
|
|
l=img.nz*(i+img.nx*j)+k;
|
|
img.z[l]=(float) raw_image[l];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Free allocated memory
|
|
free(row_pointer[0]);
|
|
free(raw_image);
|
|
|
|
// Close file
|
|
fclose(file);
|
|
|
|
return img;
|
|
}
|
|
|
|
// Write jpg
|
|
void write_jpg(char *filename,struct jpeg_image img)
|
|
{
|
|
int i,j,k,l,m;
|
|
struct jpeg_compress_struct cinfo;
|
|
struct jpeg_error_mgr jerr;
|
|
JSAMPROW row_pointer[1];
|
|
FILE *outfile;
|
|
unsigned char *raw_image=NULL;
|
|
|
|
outfile=fopen(filename,"wb");
|
|
cinfo.err=jpeg_std_error(&jerr);
|
|
jpeg_create_compress(&cinfo);
|
|
jpeg_stdio_dest(&cinfo,outfile);
|
|
cinfo.image_width=img.nx;
|
|
cinfo.image_height=img.ny;
|
|
cinfo.input_components=3;
|
|
cinfo.in_color_space=JCS_RGB;
|
|
jpeg_set_defaults(&cinfo);
|
|
jpeg_start_compress(&cinfo,TRUE);
|
|
|
|
// Allocate memory
|
|
raw_image=(unsigned char *) malloc(cinfo.image_width*cinfo.image_height*cinfo.input_components);
|
|
|
|
// Fill image
|
|
for (i=0;i<img.nx;i++) {
|
|
for (j=0;j<img.ny;j++) {
|
|
for (k=0;k<img.nz;k++) {
|
|
l=img.nz*(i+img.nx*j)+k;
|
|
raw_image[l]=(unsigned char) img.z[l];
|
|
}
|
|
}
|
|
}
|
|
|
|
while(cinfo.next_scanline<cinfo.image_height) {
|
|
row_pointer[0]=&raw_image[cinfo.next_scanline*cinfo.image_width*cinfo.input_components];
|
|
jpeg_write_scanlines(&cinfo,row_pointer,1);
|
|
}
|
|
jpeg_finish_compress(&cinfo);
|
|
jpeg_destroy_compress(&cinfo);
|
|
fclose(outfile);
|
|
|
|
return;
|
|
}
|
|
|