xref: /dports/astro/astrometry/astrometry.net-0.85/solver/codefile.c

/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

#include <string.h>
#include <unistd.h>
#include <assert.h>

#include "codefile.h"
#include "starutil.h"
#include "ioutils.h"
#include "fitsioutils.h"
#include "errors.h"
#include "quad-utils.h"

void quad_write(codefile_t* codes, quadfile_t* quads,
                unsigned int* quad, startree_t* starkd,
                int dimquads, int dimcodes) {
    double code[DCMAX];
    quad_compute_code(quad, dimquads, starkd, code);
    quad_enforce_invariants(quad, code, dimquads, dimcodes);
    codefile_write_code(codes, code);
    quadfile_write_quad(quads, quad);
}

void quad_write_const(codefile_t* codes, quadfile_t* quads,
                      const unsigned int* quad, startree_t* starkd,
                      int dimquads, int dimcodes) {
    int k;
    unsigned int quadcopy[DQMAX];
    for (k=0; k<dimquads; k++)
        quadcopy[k] = quad[k];
    quad_write(codes, quads, quadcopy, starkd, dimquads, dimcodes);
}

void codefile_compute_field_code(const double* xy, double* code, int dimquads) {
    double Ax, Ay;
    double Bx, By;
    double ABx, ABy;
    double scale, invscale;
    double costheta, sintheta;
    int i;

    Ax = xy[2*0 + 0];
    Ay = xy[2*0 + 1];
    Bx = xy[2*1 + 0];
    By = xy[2*1 + 1];
    ABx = Bx - Ax;
    ABy = By - Ay;
    scale = (ABx * ABx) + (ABy * ABy);
    invscale = 1.0 / scale;
    costheta = (ABy + ABx) * invscale;
    sintheta = (ABy - ABx) * invscale;

    // starting with star C...
    for (i=2; i<dimquads; i++) {
        double Cx, Cy;
        double x, y;
        Cx = xy[2*i + 0];
        Cy = xy[2*i + 1];
        Cx -= Ax;
        Cy -= Ay;
        x =  Cx * costheta + Cy * sintheta;
        y = -Cx * sintheta + Cy * costheta;
        code[2*(i-2)+0] = x;
        code[2*(i-2)+1] = y;
    }
}

void codefile_compute_star_code(const double* starxyz, double* code, int dimquads) {
    quad_compute_star_code(starxyz, code, dimquads);
}


#define CHUNK_CODES 0

static fitsbin_chunk_t* codes_chunk(codefile_t* cf) {
    return fitsbin_get_chunk(cf->fb, CHUNK_CODES);
}

static int callback_read_header(fitsbin_t* fb, fitsbin_chunk_t* chunk) {
    qfits_header* primheader = fitsbin_get_primary_header(fb);
    codefile_t* cf = chunk->userdata;

    cf->dimcodes = qfits_header_getint(primheader, "DIMCODES", 4);
    cf->numcodes = qfits_header_getint(primheader, "NCODES", -1);
    cf->numstars = qfits_header_getint(primheader, "NSTARS", -1);
    cf->index_scale_upper = qfits_header_getdouble(primheader, "SCALE_U", -1.0);
    cf->index_scale_lower = qfits_header_getdouble(primheader, "SCALE_L", -1.0);
    cf->indexid = qfits_header_getint(primheader, "INDEXID", 0);
    cf->healpix = qfits_header_getint(primheader, "HEALPIX", -1);
    cf->hpnside = qfits_header_getint(primheader, "HPNSIDE", 1);

    if ((cf->numcodes == -1) || (cf->numstars == -1) ||
        (cf->index_scale_upper == -1.0) || (cf->index_scale_lower == -1.0)) {
        ERROR("Couldn't find NCODES or NSTARS or SCALE_U or SCALE_L entries in FITS header");
        return -1;
    }
    if (fits_check_endian(primheader)) {
        ERROR("File was written with the wrong endianness");
        return -1;
    }
    chunk->itemsize = cf->dimcodes * sizeof(double);
    chunk->nrows = cf->numcodes;
    return 0;
}

static codefile_t* new_codefile(const char* fn, anbool writing, anbool inmem) {
    fitsbin_chunk_t chunk;
    codefile_t* cf = calloc(1, sizeof(codefile_t));
    if (!cf) {
        SYSERROR("Couldn't calloc a codefile struct");
        return NULL;
    }
    cf->healpix = -1;
    cf->hpnside = 1;

    if (inmem) {
        cf->fb = fitsbin_open_in_memory();
    } else {
        if (writing)
            cf->fb = fitsbin_open_for_writing(fn);
        else
            cf->fb = fitsbin_open(fn);
    }
    if (!cf->fb) {
        ERROR("Failed to create fitsbin");
        return NULL;
    }

    fitsbin_chunk_init(&chunk);
    chunk.tablename = "codes";
    chunk.required = 1;
    chunk.callback_read_header = callback_read_header;
    chunk.userdata = cf;
    fitsbin_add_chunk(cf->fb, &chunk);
    fitsbin_chunk_clean(&chunk);

    return cf;
}

void codefile_get_code(const codefile_t* cf, int codeid, double* code) {
    int i;
    if (codeid >= cf->numcodes) {
        ERROR("Requested codeid %i, but number of codes is %i", codeid, cf->numcodes);
        assert(codeid < cf->numcodes);
        // just carry on - we'll probably segfault.
    }
    for (i=0; i<cf->dimcodes; i++)
        code[i] = cf->codearray[codeid * cf->dimcodes + i];
}

int codefile_close(codefile_t* cf) {
    int rtn;
    if (!cf) return 0;
    rtn = fitsbin_close(cf->fb);
    free(cf);
    return rtn;
}

codefile_t* codefile_open(const char* fn) {
    codefile_t* cf = NULL;

    cf = new_codefile(fn, FALSE, FALSE);
    if (!cf)
        goto bailout;
    if (fitsbin_read(cf->fb)) {
        ERROR("Failed to open codes file");
        goto bailout;
    }
    cf->codearray = codes_chunk(cf)->data;
    return cf;

 bailout:
    if (cf)
        codefile_close(cf);
    return NULL;
}

static codefile_t* open_for_writing(const char* fn) {
    codefile_t* cf;
    qfits_header* hdr;
    if (fn)
        cf = new_codefile(fn, TRUE, FALSE);
    else
        cf = new_codefile(fn, TRUE, TRUE);
    if (!cf)
        goto bailout;
    // default
    cf->dimcodes = 4;

    // add default values to header
    hdr = codefile_get_header(cf);
    fits_add_endian(hdr);
    qfits_header_add(hdr, "AN_FILE", "CODE", "This file lists the code for each quad.", NULL);
    qfits_header_add(hdr, "NCODES", "0", "", NULL);
    qfits_header_add(hdr, "NSTARS", "0", "", NULL);
    fits_header_add_int(hdr, "DIMCODES", cf->dimcodes, "");
    qfits_header_add(hdr, "SCALE_U", "0.0", "", NULL);
    qfits_header_add(hdr, "SCALE_L", "0.0", "", NULL);
    qfits_header_add(hdr, "INDEXID", "0", "", NULL);
    qfits_header_add(hdr, "HEALPIX", "-1", "", NULL);
    qfits_header_add(hdr, "HPNSIDE", "1", "", NULL);
    fits_add_long_comment(hdr, "The first extension contains the codes "
                          "stored as %i native-endian doubles.  "
                          "(the quad location in %i-D code space)", cf->dimcodes, cf->dimcodes);
    return cf;
 bailout:
    if (cf)
        codefile_close(cf);
    return NULL;
}

codefile_t* codefile_open_for_writing(const char* fn) {
    if (!fn) {
        ERROR("Non-NULL filename required");
        return NULL;
    }
    return open_for_writing(fn);
}

codefile_t* codefile_open_in_memory() {
    return open_for_writing(NULL);
}

int codefile_switch_to_reading(codefile_t* cf) {
    if (codefile_fix_header(cf)) {
        ERROR("Failed to fix codes header");
        return -1;
    }
    if (fitsbin_switch_to_reading(cf->fb)) {
        ERROR("Failed to switch to read mode");
        return -1;
    }
    if (fitsbin_read(cf->fb)) {
        ERROR("Failed to open codes file");
        return -1;
    }
    cf->codearray = codes_chunk(cf)->data;
    return 0;
}

int codefile_write_header(codefile_t* cf) {
    fitsbin_t* fb = cf->fb;
    fitsbin_chunk_t* chunk = codes_chunk(cf);
    chunk->itemsize = cf->dimcodes * sizeof(double);
    chunk->nrows = cf->numcodes;

    if (fitsbin_write_primary_header(fb) ||
        fitsbin_write_chunk_header(fb, chunk)) {
        ERROR("Failed to write codefile header");
        return -1;
    }
    return 0;
}

int codefile_fix_header(codefile_t* cf) {
    qfits_header* hdr;
    fitsbin_t* fb = cf->fb;
    fitsbin_chunk_t* chunk = codes_chunk(cf);
    chunk->itemsize = cf->dimcodes * sizeof(double);
    chunk->nrows = cf->numcodes;

    hdr = codefile_get_header(cf);

    // fill in the real values...
    fits_header_mod_int(hdr, "DIMCODES", cf->dimcodes, "Number of values in a code.");
    fits_header_mod_int(hdr, "NCODES", cf->numcodes, "Number of codes.");
    fits_header_mod_int(hdr, "NSTARS", cf->numstars, "Number of stars.");
    fits_header_mod_double(hdr, "SCALE_U", cf->index_scale_upper, "Upper-bound index scale (radians).");
    fits_header_mod_double(hdr, "SCALE_L", cf->index_scale_lower, "Lower-bound index scale (radians).");
    fits_header_mod_int(hdr, "INDEXID", cf->indexid, "Index unique ID.");
    fits_header_mod_int(hdr, "HEALPIX", cf->healpix, "Healpix of this index.");
    fits_header_mod_int(hdr, "HPNSIDE", cf->hpnside, "Nside of the healpixelization");

    if (fitsbin_fix_primary_header(fb) ||
        fitsbin_fix_chunk_header(fb, chunk)) {
        ERROR("Failed to fix code header");
        return -1;
    }
    return 0;
}

int codefile_write_code(codefile_t* cf, double* code) {
    fitsbin_chunk_t* chunk = codes_chunk(cf);
    if (fitsbin_write_item(cf->fb, chunk, code)) {
        ERROR("Failed to write code");
        return -1;
    }
    cf->numcodes++;
    return 0;
}

qfits_header* codefile_get_header(const codefile_t* cf) {
    return fitsbin_get_primary_header(cf->fb);
}

int codefile_dimcodes(const codefile_t* cf) {
    return cf->dimcodes;
}