#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

# include "psphot.h"

int psphotSaveImage (psMetadata *header, psImage *image, char *filename) {

    psFits *fits = psFitsOpen (filename, "w");
    psFitsWriteImage (fits, NULL, image, 0, NULL);
    psFitsClose (fits);
    return (TRUE);
}

bool psphotDumpMoments (psMetadata *recipe, psArray *sources) {

    bool status;

    // optional dump of all rough source data
    char *output = psMetadataLookupStr (&status, recipe, "MOMENTS_OUTPUT_FILE");
    if (!status) return false;
    if (output == NULL) return false;
    if (output[0] == 0) return false;

    pmMomentsWriteText (sources, output);
    return true;
}

bool psphotDumpSource (pmSource *source, char *name) {

    FILE *f = fopen (name, "w");
    if (f == NULL) psAbort ("psphot", "can't open file");

    for (int i = 0; i < source->pixels->numRows; i++) {
        for (int j = 0; j < source->pixels->numCols; j++) {
            // skip masked points
            if (source->mask->data.U8[i][j]) {
                continue;
            }
            // skip zero-weight points
            if (source->weight->data.F32[i][j] == 0) {
                continue;
            }

	    fprintf (f, "%d %d %f %f %d\n", 
		     (j + source->pixels->col0),
		     (i + source->pixels->row0),
		     source->pixels->data.F32[i][j],
		     1.0 / source->weight->data.F32[i][j],
		     source->mask->data.U8[i][j]);
        }
    }
    fclose (f);
    return true;
}

// these values are saved in an output header stub - they are added to either the
// PHU header (CMP) or the MEF table header (CMF)
psMetadata *psphotDefineHeader (psMetadata *recipe) {

    psMetadata *header = psMetadataAlloc ();

    // write necessary information to output header
    psMetadataItemTransfer (header, recipe, "ZERO_PT");
    psMetadataItemTransfer (header, recipe, "PHOTCODE");

    psMetadataItemTransfer (header, recipe, "APMIFIT");
    psMetadataItemTransfer (header, recipe, "DAPMIFIT");
    psMetadataItemTransfer (header, recipe, "NAPMIFIT");
    psMetadataItemTransfer (header, recipe, "SKYBIAS");
    psMetadataItemTransfer (header, recipe, "SKYSAT");
    
    // PSF model parameters (shape values for image center)
    psMetadataItemTransfer (header, recipe, "NPSFSTAR");
    psMetadataItemTransfer (header, recipe, "APLOSS");
    psMetadataItemTransfer (header, recipe, "FWHM_X");
    psMetadataItemTransfer (header, recipe, "FWHM_Y");
    psMetadataItemTransfer (header, recipe, "ANGLE");

    // XXX these need to be defined from elsewhere
    psMetadataAdd (header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
    psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);

    // sky background model statistics
    psMetadataItemTransfer (header, recipe, "SKY_MEAN");
    psMetadataItemTransfer (header, recipe, "SKY_STDV");
    psMetadataItemTransfer (header, recipe, "SKY_MIN");
    psMetadataItemTransfer (header, recipe, "SKY_MAX");
    psMetadataItemTransfer (header, recipe, "SKY_NX");
    psMetadataItemTransfer (header, recipe, "SKY_NY");

    return header;
}

# if (0)
// output functions: we have several fixed modes (sx, obj, cmp)
void psphotOutput (pmReadout *readout, psMetadata *arguments) {

    bool status;
    char *outputFile = NULL;

    psMetadata *header = pmReadoutGetHeader (readout);

    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    // sample PSF images??
    if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);
    if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);

    if (psfFile != NULL) {
	psMetadata *psfData = pmPSFtoMD (NULL, psf);
	psMetadataConfigWrite (psfData, psfFile);
	psFree (psfData);
    }
}

psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) {

    psRegion region = {x - dx, x + dx, y - dy, y + dy};
    psImage *view = psImageSubset (image, region);
    psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
    psImageInit (sample, 0);
    modelEXT->params->data.F32[2] = x;
    modelEXT->params->data.F32[3] = y;
    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    pmModelAdd (sample, NULL, modelPSF, false, false);
    psFree (modelPSF);
    return (sample);
}

bool psphotSamplePSFs (pmPSF *psf, psImage *image, char *output) {

    // make sample PSFs for 4 corners and the center
    psImage *sample;

    // optional dump of all rough source data
    if (output[0] == 0) return false;

    pmModel *modelEXT = pmModelAlloc (psf->type);
    modelEXT->params->data.F32[0] = 0;
    modelEXT->params->data.F32[1] = 1;

    psFits *fits = psFitsOpen (output, "w");

    // the centers are in parent coordinates; they do not need to correspond to valid pixels...
    sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25);
    psFitsWriteImage (fits, NULL, sample, 0);
    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
    psFitsWriteImage (fits, NULL, sample, 0);
    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25);
    psFitsWriteImage (fits, NULL, sample, 0);
    sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25);
    psFitsWriteImage (fits, NULL, sample, 0);
    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
    psFitsWriteImage (fits, NULL, sample, 0);

    psFitsClose (fits);
    return (TRUE);
}
# endif
