356 lines
8.0 KiB
C
356 lines
8.0 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;
|
|
}
|