/* $Id: create_mosaic.c,v 1.10 2002/01/15 10:22:17 fors Exp $ * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * COPYRIGHT (c) 2001 European Southern Observatory * LICENSE: GNU General Public License version 2 or later * * PROJECT: VLT Data Flow System * AUTHOR: Ralf Palsa -- ESO/DMD/DPG * SUBSYSTEM: Instrument pipelines * * PURPOSE: * DESCRIPTION: * * $Name: fsmosaic-1_0 $ * $Revision: 1.10 $ * ---------------------------------------------------------------------------- */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #ifdef HAVE_DIRENT_H # include #else # define direct dirent # ifdef HAVE_SYS_NDIR_H # include # endif # ifdef HAVE_SYS_DIR_H # include # endif # ifdef HAVE_NDIR_H # include # endif #endif #include #include "image.h" #include "qfits_ext.h" #include "overscan.h" #include "align.h" #include "create_mosaic.h" #include "file_utils.h" /* * Verbosity flag */ int verbose = 0; /* * Structure defining selection criteria for file selection. */ struct selection { char instrument[FITS_CARD_SIZE + 1]; char chip_id[FITS_CARD_SIZE + 1]; char product_id[FITS_CARD_SIZE + 1]; int is_product; int n_chips; double mjd; }; /* * @memo * Search for the appropriate master or slave file. * * @return The function returns the name of the master or slave file, if no * error occurred, otherwise a #NULL# pointer is returned. * * @param scandir Directory to scan * @param query Selection criteria * * @doc * The function scans the directory given by \textbf{scandir} for all * FITS files. For each FITS file found, the header keywords providing * the instrument name, the chip identification string, the number of * chips and the Modified Julian date match the ones provided by the * structure members \textbf{instrument}, \textbf{chip\_id}, * \textbf{n\_chips}, and \textbf{mjd} of \textbf{query}. * * In addition to the above comparision the value of the structure member * \textbf{product_id} is compared to the header keyword providing the * pipeline product tag, if the structure member \textbf{is\_product} is * different from 0. * * The function returns the name of the first FITS file for which all of * the above comparisons succeed. * * The modified julian date is considered to be the same if the difference * between the member \textbf{mjd} and the modified julian date of the * current FITS file is less than #TINY#. * * @author R. Palsa */ static char *select_file(const char *scandir, struct selection *query) { static char path[PATHNAME_MAX + 1]; char *name = 0; size_t name_max, path_max; int found; int valid_instrument, valid_chip_id, valid_n_chips, valid_mjd; int valid_product; struct dirent *entry; DIR *dir; /* * Get maximum possible length of a file name in the search directory, * as defined by the system. */ #ifdef HAVE_PATHCONF name_max = pathconf(scandir, _PC_NAME_MAX); #else # ifdef NAME_MAX name_max = NAME_MAX; # else name_max = PATHNAME_MAX; # endif #endif if ((path_max = strlen(scandir)) > PATHNAME_MAX - name_max - 1) return 0; if (!(dir = opendir(scandir))) return 0; /* * Build file name root. */ strncpy(path, scandir, path_max); path[path_max] = '/'; name = path + path_max + 1; /* * Scan source directory for FITS files and compare the header entries * specifying the instrument name, modified julian date and CCD * identification string. */ found = 0; entry = 0; while (!found && (entry = readdir(dir))) { strncpy(name, entry->d_name, name_max + 1); if (is_fits_file(path)) { char *s; valid_instrument = 0; valid_chip_id = 0; valid_n_chips = 0; valid_mjd = 0; valid_product = 1; if ((s = qfits_query_hdr(path, INSTRUMENT)) && !strcmp(qfits_pretty_string(s), query->instrument)) valid_instrument = 1; if ((s = qfits_query_hdr(path, CHIP1_ID)) && !strcmp(qfits_pretty_string(s), query->chip_id)) valid_chip_id = 1; if ((s = qfits_query_hdr(path, N_CHIPS))) if (atoi(s) == query->n_chips) valid_n_chips = 1; if ((s = qfits_query_hdr(path, MJD))) { double mjd = atof(s); if (fabs(mjd - query->mjd) <= TINY) valid_mjd = 1; } /* * Honor product tag only in case of a pipeline processed frame. */ if (query->is_product) { if (!(s = qfits_query_hdr(path, PRO_CATEGORY)) || strcmp(qfits_pretty_string(s), query->product_id)) valid_product = 0; } if (valid_instrument && valid_chip_id && valid_n_chips && valid_mjd && valid_product) found = 1; else found = 0; } } closedir(dir); if (!found || !entry) return 0; return path; } /* * @memo * Merge the master and slave FITS headers. * * @return The function returns the pointer to the merged header, if no error * occurred, otherwise the return value is #NULL#. * * @param header1 First header to be merged. * @param header2 Second header to be merged. * @param naxis1 Size of merged frame in X-direction. * @param naxis2 Size of merged frame in Y-direction. * @param xoffset Offset of the reference pixel along X. * @param yoffset Offset of the reference pixel along Y. * * @doc * The function merges the input FITS headers \textbf{header1} and * \textbf{header2} into a new output header, correcting the reference * pixel location. The reference pixel of the merged header is the * location of the reference pixel of \textbf{header1} with the offsets * along the X and Y axis added given by \textbf{xoffset} and * \textbf{yoffset} respectively. * * The merged header is basically \textbf{header1} updated with the * renamed, detector related keywords from \textbf{header2}. * * @author R. Palsa */ static qfits_header *merge_headers(qfits_header *header1, qfits_header *header2, int naxis1, int naxis2, double xoffset, double yoffset) { char card[FITS_CARD_SIZE + 1]; char *value, *comment; unsigned int outputs1, outputs2; double crpix1, crpix2; qfits_header *merged_header; /* * Initialize merged header with header1. */ if (!(merged_header = qfits_header_copy(header1))) return 0; /* * Correct reference pixel location in header1 to account for the * new image size. */ /* FIXME: * This assumes that the reference pixel is located on the CCD which * need not be the case. Reference pixels laying outside the CCD are * allowed, i.e. negative values of CRPIXi may identify valid locations. */ crpix1 = qfits_header_getdouble(header1, "CRPIX1", -1); crpix2 = qfits_header_getdouble(header1, "CRPIX2", -1); if (crpix1 < 0 || crpix2 < 0) return 0; crpix1 += xoffset; crpix2 += yoffset; sprintf(card, "%.1f", crpix1); fits_header_set_value(merged_header, "CRPIX1", card); sprintf(card, "%.1f", crpix2); fits_header_set_value(merged_header, "CRPIX2", card); /* * Copy detector keywords from header2 to header1 renaming the * keywords according to the DICB rules. */ /* * CHIP keywords */ fits_header_set_value(merged_header, N_CHIPS, "2"); fits_header_set_value(merged_header, CHIP1_INDEX, "1"); fits_header_set_value(merged_header, CHIP1_X, "1"); fits_header_set_value(merged_header, CHIP1_Y, "2"); comment = qfits_header_getcom(header2, CHIP1_ID); qfits_header_add(merged_header, CHIP2_ID, CHIP_SLAVE, comment, 0); value = qfits_header_getstr(header2, CHIP1_NAME); comment = qfits_header_getcom(header2, CHIP1_NAME); qfits_header_add(merged_header, CHIP2_NAME, value, comment, 0); value = qfits_header_getstr(header2, CHIP1_DATE); comment = qfits_header_getcom(header2, CHIP1_DATE); qfits_header_add(merged_header, CHIP2_DATE, value, comment, 0); comment = qfits_header_getcom(header2, CHIP1_INDEX); qfits_header_add(merged_header, CHIP2_INDEX, "2", comment, 0); comment = qfits_header_getcom(header2, CHIP1_X); qfits_header_add(merged_header, CHIP2_X, "1", comment, 0); comment = qfits_header_getcom(header2, CHIP1_Y); qfits_header_add(merged_header, CHIP2_Y, "1", comment, 0); sprintf(card, "%d", qfits_header_getint(header2, CHIP1_NX, -1)); comment = qfits_header_getcom(header2, CHIP1_NX); qfits_header_add(merged_header, CHIP2_NX, card, comment, 0); sprintf(card, "%d", qfits_header_getint(header2, CHIP1_NY, -1)); comment = qfits_header_getcom(header2, CHIP1_NY); qfits_header_add(merged_header, CHIP2_NY, card, comment, 0); sprintf(card, "%.1f", qfits_header_getdouble(header2, CHIP1_PSZX, 0)); comment = qfits_header_getcom(header2, CHIP1_PSZX); qfits_header_add(merged_header, CHIP2_PSZX, card, comment, 0); sprintf(card, "%.1f", qfits_header_getdouble(header2, CHIP1_PSZY, 0)); comment = qfits_header_getcom(header2, CHIP1_PSZY); qfits_header_add(merged_header, CHIP2_PSZY, card, comment, 0); /* * OUTPUT keywords */ if ((outputs1 = qfits_header_getint(header1, N_OUTPUTS, 0)) < 1) { qfits_header_destroy(merged_header); return 0; } if ((outputs2 = qfits_header_getint(header2, N_OUTPUTS, 0)) < 1) { qfits_header_destroy(merged_header); return 0; } sprintf(card, "%d", outputs1 + outputs2); fits_header_set_value(merged_header, N_OUTPUTS, card); fits_header_set_value(merged_header, OUTPUT1_INDEX, "1"); fits_header_set_value(merged_header, OUTPUT1_CHIP, "1"); qfits_header_del(merged_header, OUTPUT1_X); qfits_header_del(merged_header, OUTPUT1_Y); qfits_header_del(merged_header, OUTPUT1_NX); qfits_header_del(merged_header, OUTPUT1_NY); fits_header_set_value(merged_header, OUTPUT1_PRSCY, "0"); sprintf(card, "%-8s", qfits_header_getstr(header2, OUTPUT1_ID)); comment = qfits_header_getcom(header2, OUTPUT1_ID); qfits_header_add(merged_header, OUTPUT2_ID, card, comment, 0); sprintf(card, "%-8s", qfits_header_getstr(header2, OUTPUT1_NAME)); comment = qfits_header_getcom(header2, OUTPUT1_NAME); qfits_header_add(merged_header, OUTPUT2_NAME, card, comment, 0); comment = qfits_header_getcom(header2, OUTPUT1_INDEX); qfits_header_add(merged_header, OUTPUT2_INDEX, "2", comment, 0); comment = qfits_header_getcom(header2, OUTPUT1_CHIP); qfits_header_add(merged_header, OUTPUT2_CHIP, "2", comment, 0); sprintf(card, "%d", qfits_header_getint(header2, OUTPUT1_PRSCX, -1)); comment = qfits_header_getcom(header2, OUTPUT1_PRSCX); qfits_header_add(merged_header, OUTPUT2_PRSCX, card, comment, 0); sprintf(card, "%d", qfits_header_getint(header2, OUTPUT1_OVSCX, -1)); comment = qfits_header_getcom(header2, OUTPUT1_OVSCX); qfits_header_add(merged_header, OUTPUT2_OVSCX, card, comment, 0); sprintf(card, "%d", qfits_header_getint(header2, OUTPUT1_PRSCY, -1)); comment = qfits_header_getcom(header2, OUTPUT1_PRSCY); qfits_header_add(merged_header, OUTPUT2_PRSCY, card, comment, 0); sprintf(card, "%d", 0); comment = qfits_header_getcom(header2, OUTPUT1_OVSCY); qfits_header_add(merged_header, OUTPUT2_OVSCY, card, comment, 0); sprintf(card, "%.2f", qfits_header_getdouble(header2, OUTPUT1_CONAD, 0)); comment = qfits_header_getcom(header2, OUTPUT1_CONAD); qfits_header_add(merged_header, OUTPUT2_CONAD, card, comment, 0); sprintf(card, "%.2f", qfits_header_getdouble(header2, OUTPUT1_GAIN, 0)); comment = qfits_header_getcom(header2, OUTPUT1_GAIN); qfits_header_add(merged_header, OUTPUT2_GAIN, card, comment, 0); sprintf(card, "%.2f", qfits_header_getdouble(header2, OUTPUT1_RON, 0)); comment = qfits_header_getcom(header2, OUTPUT1_RON); qfits_header_add(merged_header, OUTPUT2_RON, card, comment, 0); /* * WINDOW keywords */ sprintf(card, "%d", naxis1); fits_header_set_value(merged_header, WINDOW1_NX, card); sprintf(card, "%d", naxis2); fits_header_set_value(merged_header, WINDOW1_NY, card); /* * Delete original DPR keywords from the merged header and write * product tag instead. */ qfits_header_del(merged_header, DPR_CATEGORY); qfits_header_del(merged_header, DPR_TYPE); qfits_header_del(merged_header, DPR_TECHNIQUE); qfits_header_add(merged_header, PRO_CATEGORY, MOSAIC_ID, "Product tag", 0); /* * Delete other keys */ qfits_header_del(merged_header, EXTNAME); qfits_header_del(merged_header, ORIGFILE); qfits_header_del(merged_header, ARCFILE); qfits_header_del(merged_header, CHECKSUM); return merged_header; } /** * @memo * Check if a file contains a merged image. * * @return The function returns 1 if the file contains a merged image. * Otherwise the return value is 0. * * @param filename Name of the file to be checked. * * @doc * The function checks that the input file given by \textbf{filename} is * a FITS file. To be recognized as a mosaic the header of the file * \textbf{filename} must contain the product tag #FORS2_MERGED_IMAGE#, * the number of chips must be 2 and the chip identifiers must be the * appropriate master and slave tags. * * @author R.Palsa */ int is_mosaic(const char *filename) { char *s; if (!is_fits_file((char *)filename)) return -1; s = qfits_query_hdr((char *)filename, PRO_CATEGORY); if (!s || strcmp(qfits_pretty_string(s), MOSAIC_ID)) return 0; s = qfits_query_hdr((char *)filename, N_CHIPS); if (!s || atoi(s) != 2) return 0; /* * Do a loose check on the chips identifiers. A particular order * of master and slave is not required. */ s = qfits_query_hdr((char *)filename, CHIP1_ID); if (!s || (strcmp(qfits_pretty_string(s), CHIP_SLAVE) && strcmp(qfits_pretty_string(s), CHIP_MASTER))) return 0; s = qfits_query_hdr((char *)filename, CHIP2_ID); if (!s || (strcmp(qfits_pretty_string(s), CHIP_SLAVE) && strcmp(qfits_pretty_string(s), CHIP_MASTER))) return 0; return 1; } /** * @memo * Create a mosaic from a FORS2 master and slave frame. * * @return The function returns #EXIT_SUCCESS# if no error occured, otherwise * one of the return codes defined in create_mosaic.h. * * @param input_file Filename of the input master or slave frame. * @param output_file Filename of the output mosaic frame. * * @doc * TBD * * @author R. Palsa */ int create_mosaic(const char *input_file, const char *output_file) { const char *fctid = "create_mosaic"; char *s, *dir; char *master_file, *slave_file; int master_flag = 0; int binx, biny; int bits_per_pixel; double xshift, yshift, angle; double xsize, ysize; double xoffset, yoffset; image_t *master_image, *slave_image, *mosaic; struct selection query; fits_info_t *info; qfits_header *master_header, *slave_header, *merged_header; /* * Get instrument name and modified julian date from * the input FITS file and setup the file selection criteria * for the directory scan. */ if (!(s = qfits_query_hdr((char *)input_file, INSTRUMENT))) { fprintf(stderr, "%s: %s: Missing FITS keyword '%s'!\n", fctid, input_file, INSTRUMENT); return EXIT_NO_INSTRUMENT; } else { s = qfits_pretty_string(s); if (strcmp(s, INSTRUMENT_ID)) { fprintf(stderr, "%s: %s: Invalid instrument '%s'!\n", fctid, input_file, s); return EXIT_UNKNOWN_INSTRUMENT; } else strncpy(query.instrument, s, FITS_CARD_SIZE); } if (!(s = qfits_query_hdr((char *)input_file, N_CHIPS))) { fprintf(stderr, "%s: %s: Missing FITS keyword '%s'!\n", fctid, input_file, N_CHIPS); return EXIT_MISSING_KEYWORD; } else { query.n_chips = atoi(s); if (query.n_chips != 1) { fprintf(stderr, "%s: %s: Invalid number of chips!\n", fctid, input_file); return EXIT_WRONG_KEY_VALUE; } } if (!(s = qfits_query_hdr((char *)input_file, MJD))) { fprintf(stderr, "%s: %s: Missing FITS keyword '%s'!\n", fctid, input_file, MJD); return EXIT_MISSING_KEYWORD; } else query.mjd = atof(s); /* * Check if the product tag keyword as created by the pipeline for * processed frames. Its value will be used to discriminate between * different products of a single recipe run. If it is not present * it is not an error, but it is assumed that the input frame is a * raw frame. */ if ((s = qfits_query_hdr((char *)input_file, PRO_CATEGORY))) { query.is_product = 1; strncpy(query.product_id, qfits_pretty_string(s), FITS_CARD_SIZE); } else query.is_product = 0; /* * Determine the master and slave frame. The current scheme assumes that * the input file passed to the function is located in the current working * directory and that the secondary file is present in this directory too. * To find the secondary file, the current directory is scanned for all * FITS files. Among these FITS files the one having the same INSTRUME * entry and the same MJD-OBS entry is selected as the secondary file. * The decision which of the two file is the master and slave respectively * is based on the CCD id's in the header. */ dir = dir_name(input_file); if (!(s = qfits_query_hdr((char *)input_file, CHIP1_ID))) { fprintf(stderr, "%s: %s: Missing FITS keyword '%s'!\n", fctid, input_file, CHIP1_ID); return EXIT_MISSING_KEYWORD; } else { char *chip_id = strdup(qfits_pretty_string(s)); if (!strcmp(chip_id, CHIP_MASTER)) { master_file = (char *)input_file; master_flag = 1; strncpy(query.chip_id, CHIP_SLAVE, FITS_CARD_SIZE); if (!(s = select_file(dir, &query))) { fprintf(stderr, "%s: %s: No appropriate slave file found for " "master!\n", fctid, master_file); free(chip_id); return EXIT_NO_SLAVE_FOUND; } else slave_file = strdup(s); } else { if (!strcmp(chip_id, CHIP_SLAVE)) { slave_file = (char *)input_file; master_flag = 0; strncpy(query.chip_id, CHIP_MASTER, FITS_CARD_SIZE); if (!(s = select_file(dir, &query))) { fprintf(stderr, "%s: %s: No appropriate master file found for " "slave!\n", fctid, slave_file); free(chip_id); return EXIT_NO_MASTER_FOUND; } else master_file = strdup(s); } else { fprintf(stderr, "%s: %s: Unknown CCD ID '%s'!\n", fctid, input_file, chip_id); free(chip_id); return EXIT_UNKNOWN_CCD_ID; } } free(chip_id); } if (verbose) { fprintf(stdout, "Master file: %s\n", master_file); fprintf(stdout, "Slave file: %s\n", slave_file); } /* * Get the original pixel format of the master file. */ if (!(info = fits_get_file_info(master_file))) { fprintf(stderr, "%s: %s: Cannot read basic FITS file information!\n", fctid, master_file); fits_info_del(info); if (master_flag) free(slave_file); else free(master_file); return EXIT_FITS_ERROR; } bits_per_pixel = info->bits_per_pixel; fits_info_del(info); /* * Load image and header data. */ if (verbose) fprintf(stdout, "Loading images ..."); if (!(master_header = qfits_header_read(master_file))) { fprintf(stderr, "%s: %s: Cannot read FITS header!\n", fctid, master_file); return EXIT_FITS_ERROR; } if (!(master_image = fits_load_image(master_file))) { fprintf(stderr, "%s: %s: Cannot load image!\n", fctid, master_file); qfits_header_destroy(master_header); return EXIT_FITS_ERROR; } if (!(slave_header = qfits_header_read(slave_file))) { fprintf(stderr, "%s: %s: Cannot read FITS header!\n", fctid, slave_file); image_del(master_image); qfits_header_destroy(master_header); return EXIT_FITS_ERROR; } if (!(slave_image = fits_load_image(slave_file))) { fprintf(stderr, "%s: %s: Cannot load image!\n", fctid, slave_file); image_del(master_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); return EXIT_FITS_ERROR; } if (master_flag) free(slave_file); else free(master_file); if (verbose) fprintf(stdout, " done\n"); /* * Get shifts and rotation angle from the slave image header. * It is assumed that the keywords exist and have reasonable values. * In case the keywords are not found in the slave header, the shifts * and the rotation angle are set to 0., i.e. the images are just put * into the output image as they are. */ xshift = qfits_header_getdouble(slave_header, XSHIFT, 0.); yshift = qfits_header_getdouble(slave_header, YSHIFT, 0.); angle = qfits_header_getdouble(slave_header, ROTATION, 0.); if ((xsize = qfits_header_getdouble(slave_header, XSIZE, -1)) < 0.) { fprintf(stderr, "%s: Invalid X pixel size!\n", fctid); image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); return EXIT_WRONG_KEY_VALUE; } if ((ysize = qfits_header_getdouble(slave_header, YSIZE, -1)) < 0.) { fprintf(stderr, "%s: Invalid Y pixel size!\n", fctid); image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); return EXIT_WRONG_KEY_VALUE; } binx = qfits_header_getint(slave_header, BINX, 1); biny = qfits_header_getint(slave_header, BINY, 1); xshift /= xsize * binx; yshift /= ysize * biny; xshift = -xshift; yshift = -yshift; if (verbose) { fprintf(stdout, "Rotating image by %.4f degrees\n", angle); fprintf(stdout, "Shifting image by %.4f, %.4f pixels\n", xshift, yshift); } /* * Align the two images assuming that the CCDs are stacked vertically. * Remove overscan areas between the two CCDs (along y) first if necessary. */ if (trim_prescan(master_header, master_image, 1) == EXIT_FAILURE) { fprintf(stderr, "%s: Cannot remove prescan area from master!\n", fctid); image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); return EXIT_FAILURE; } if (trim_overscan(slave_header, slave_image, 1) == EXIT_FAILURE) { fprintf(stderr, "%s: Cannot remove overscan area from slave!\n", fctid); image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); return EXIT_FAILURE; } yshift -= master_image->vsize; mosaic = align_images(master_image, slave_image, angle, xshift, yshift, &xoffset, &yoffset); if (!mosaic) { image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); if (verbose) fprintf(stdout, "failed\n"); return EXIT_ALIGN_ERROR; } /* * Update the coordinate system origin in the FITS header. */ if (verbose) fprintf(stdout, "Updating FITS header ..."); if (!(merged_header = merge_headers(master_header, slave_header, mosaic->hsize, mosaic->vsize, xoffset, yoffset))) { image_del(master_image); image_del(slave_image); qfits_header_destroy(master_header); qfits_header_destroy(slave_header); if (verbose) fprintf(stdout, "failed\n"); return EXIT_FAILURE; } qfits_header_destroy(master_header); qfits_header_destroy(slave_header); if (verbose) fprintf(stdout, " done\n"); /* * Write the created mosaic to disk. */ if (verbose) fprintf(stdout, "Saving result to %s ...", output_file); if (fits_save_image(output_file, merged_header, mosaic, bits_per_pixel) == EXIT_FAILURE) { fprintf(stderr, "%s: Cannot write file '%s'\n", fctid, output_file); image_del(master_image); image_del(slave_image); image_del(mosaic); qfits_header_destroy(merged_header); return EXIT_SAVE_ERROR; } if (verbose) fprintf(stdout, " done\n"); /* * Cleanup */ image_del(master_image); image_del(slave_image); image_del(mosaic); qfits_header_destroy(merged_header); return EXIT_SUCCESS; }