Index: trunk/ppSkycell/src/Makefile.am =================================================================== --- trunk/ppSkycell/src/Makefile.am (revision 23982) +++ trunk/ppSkycell/src/Makefile.am (revision 23992) @@ -31,6 +31,7 @@ ppSkycell.c \ ppSkycellArguments.c \ + ppSkycellCamera.c \ + ppSkycellData.c \ ppSkycellLoop.c \ - ppSkycellReadout.c \ ppSkycellVersion.c Index: trunk/ppSkycell/src/ppSkycell.c =================================================================== --- trunk/ppSkycell/src/ppSkycell.c (revision 23982) +++ trunk/ppSkycell/src/ppSkycell.c (revision 23992) @@ -5,4 +5,7 @@ #include #include +#include + +#include "ppSkycell.h" int main(int argc, char *argv[]) @@ -11,5 +14,5 @@ if (!data) { psErrorStackPrint(stderr, "Unable to initialise."); - exit(PS_EXIT_CONFIG_ERROR); + return PS_EXIT_CONFIG_ERROR; } @@ -17,7 +20,22 @@ psErrorStackPrint(stderr, "Unable to parse arguments."); psFree(data); - exit(PS_EXIT_CONFIG_ERROR); + return PS_EXIT_CONFIG_ERROR; } + if (!ppSkycellCamera(data)) { + psErrorStackPrint(stderr, "Unable to parse camera configuration."); + psFree(data); + return PS_EXIT_CONFIG_ERROR; + } + if (!ppSkycellLoop(data)) { + psErrorStackPrint(stderr, "Unable to process data."); + psFree(data); + return PS_EXIT_DATA_ERROR; + } + psFree(data); + + return PS_EXIT_SUCCESS; +} + Index: trunk/ppSkycell/src/ppSkycell.h =================================================================== --- trunk/ppSkycell/src/ppSkycell.h (revision 23982) +++ trunk/ppSkycell/src/ppSkycell.h (revision 23992) @@ -2,10 +2,37 @@ #define PP_SKYCELL_H +#include +#include + +#define PPSKYCELL_RECIPE "PPSKYCELL" // Recipe name + +// Data for processing typedef struct { - psString inName; // Input filename - psString outName; // Output name + psString imagesName; // Filename with images + psString masksName; // Filename with masks + psString outRoot; // Output root name int numInputs; // Number of inputs + psImageMaskType maskVal; // Value to mask int bin1, bin2; // Binning factors pmConfig *config; // Configuration } ppSkycellData; +/// Initialise data for processing +ppSkycellData *ppSkycellDataInit(int *argc, char *argv[] // Command-line arguments + ); + +/// Parse command-line arguments +bool ppSkycellArguments(ppSkycellData *data, // Data for processing + int argc, char *argv[] // Command-line arguments + ); + +/// Parse camera configurations +bool ppSkycellCamera(ppSkycellData *data // Data for processing + ); + +/// Loop over input data, processing +bool ppSkycellLoop(ppSkycellData *data // Data for processing + ); + + +#endif Index: trunk/ppSkycell/src/ppSkycellArguments.c =================================================================== --- trunk/ppSkycell/src/ppSkycellArguments.c (revision 23982) +++ trunk/ppSkycell/src/ppSkycellArguments.c (revision 23992) @@ -29,5 +29,5 @@ { fprintf(stderr, "\nPan-STARRS skycell JPEGifier\n\n"); - fprintf(stderr, "Usage: %s INPUT.list OUTPUT_ROOT\n\n", + fprintf(stderr, "Usage: %s -images INPUT.list [-masks MASK.list] OUTPUT_ROOT\n\n", program); fprintf(stderr, "\n"); @@ -45,16 +45,23 @@ bool ppSkycellArguments(ppSkycellData *data, int argc, char *argv[]) { - assert(config); + assert(data); + assert(data->config); psMetadata *arguments = psMetadataAlloc(); // Command-line arguments - if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { + psMetadataAddStr(arguments, PS_LIST_TAIL, "-images", 0, "Filename with input images", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-masks", 0, "Filename with input masks", NULL); + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { usage(argv[0], arguments, data); } - data->inName = psStringCopy(argv[1]); - data->outRoot = psStringCopy(argv[2]); + data->imagesName = psMetadataLookupStr(NULL, arguments, "-images"); + data->masksName = psMetadataLookupStr(NULL, arguments, "-masks"); + data->outRoot = psStringCopy(argv[1]); psTrace("ppSkycell", 1, "Done reading command-line arguments\n"); psFree(arguments); + + PS_ASSERT_STRING_NON_EMPTY(data->imagesName, false); + PS_ASSERT_STRING_NON_EMPTY(data->outRoot, false); return true; Index: trunk/ppSkycell/src/ppSkycellCamera.c =================================================================== --- trunk/ppSkycell/src/ppSkycellCamera.c (revision 23982) +++ trunk/ppSkycell/src/ppSkycellCamera.c (revision 23992) @@ -9,13 +9,35 @@ #include "ppSkycell.h" +/// Read list of images +static psArray *fileList(const char *filename // Filename + ) +{ + psString input = psSlurpFilename(filename); + if (!input) { + psError(PS_ERR_IO, false, "Unable to read %s", filename); + return false; + } + psArray *inputs = psStringSplitArray(input, "\n", false); // Input filenames + psFree(input); + if (!inputs || inputs->n == 0) { + psError(PS_ERR_IO, false, "Unable to read filenames from %s", filename); + psFree(inputs); + return NULL; + } + return inputs; +} + /// Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc -static void fileList(const char *file, // The symbolic name for the file - const char *name, // The name of the file - const char *comment, // Description of the file - pmConfig *config // Configuration +static void fileArguments(const char *file, // The symbolic name for the file + const char *name, // The name of the file + const char *comment, // Description of the file + pmConfig *config // Configuration ) { psArray *files = psArrayAlloc(1); // Array with file names files->data[0] = psStringCopy(name); + if (psMetadataLookup(config->arguments, file)) { + psMetadataRemoveKey(config->arguments, file); + } psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); psFree(files); @@ -27,44 +49,55 @@ ) { - psString input = psSlurpFilename(data->inName); - if (!input) { - psError(PS_ERR_IO, false, "Unable to read %s", data->inName); + psArray *images = fileList(data->imagesName); // Image names + if (!images) { + psError(psErrorCodeLast(), false, "No images provided."); return false; } + data->numInputs = images->n; - psArray *inputs = psStringSplitArray(input, "\n", false); // Input filenames - psFree(input); - if (!inputs || inputs->n == 0) { - psError(PS_ERR_IO, false, "Unable to read filenames from %s", data->inName); - psFree(inputs); - return false; - } - - data->numInputs = inputs->n; - psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Input files", inputs); - psFree(inputs); - psMetadataAddString(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot); - - pmFPAfile *inFile = NULL; // Representative input file - for (int i = 0; i < data->numInputs; i++) { - bool found = false; // Found file definition? - pmFPAfile *file = pmFPAfileDefineSingleFromArgs(&found, data->config, "PPSKYCELL.INPUT", - "INPUTS", i); // Input file - if (!file || !found) { - psError(psErrorCodeLast(), false, "Unable to define input %d\n", i); + psArray *masks = NULL; // Mask names + if (data->masksName) { + masks = fileList(data->masksName); + if (!masks) { + psError(psErrorCodeLast(), false, "No masks provided."); + psFree(images); return false; } - if (!inFile) { - inFile = file; + if (masks->n != data->numInputs) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of images (%ld) and masks (%ld) do not match", + images->n, masks->n); + psFree(images); + psFree(masks); + return false; } } - pmFPAfile *outFile = pmFPAfileDefineFromFile(data->config, inFile, data->bin1, data->bin1, - "PPSKYCELL.JPEG1"); // Output file - if (!outFile) { + psMetadataAddStr(data->config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot); + + for (int i = 0; i < data->numInputs; i++) { + bool status = false; // Status of file definition + fileArguments("IMAGE", images->data[i], "Name of the image", data->config); + pmFPAfile *image = pmFPAfileDefineFromArgs(&status, data->config, "PPSKYCELL.IMAGE", "IMAGE"); // File + if (!status || !image) { + psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.IMAGE"); + // XXX Cleanup + return false; + } + + if (data->masksName) { + fileArguments("MASK", masks->data[i], "Name of the mask", data->config); + if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.MASK", "MASK") || !status) { + psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.MASK"); + // XXX Cleanup + return false; + } + } + } + + if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.JPEG1")) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; } - if (!pmFPAfileDefineFromFile(data->config, outFile, data->bin2, data->bin2, "PPSKYCELL.JPEG2")) { + if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.JPEG2")) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; Index: trunk/ppSkycell/src/ppSkycellLoop.c =================================================================== --- trunk/ppSkycell/src/ppSkycellLoop.c (revision 23982) +++ trunk/ppSkycell/src/ppSkycellLoop.c (revision 23992) @@ -11,13 +11,28 @@ #define BUFFER 16 // Size of buffer for projections -static psRegion *skycellRegion(psRegion *region, // Current region of skycell, or NULL - pmAstromWCS *wcs, // World Coordinate System +// List of input files +static const char *inFiles[] = { "PPSKYCELL.IMAGE", "PPSKYCELL.MASK", NULL }; + +// List of output files +//static const char *outFiles[] = { "PPSKYCELL.JPEG1", "PPSKYCELL.JPEG2", NULL }; + + +static void regionMinMax(psRegion *base,// Base region; modified + const psRegion *compare // Comparison region + ) +{ + base->x0 = PS_MIN(base->x0, compare->x0); + base->x1 = PS_MAX(base->x1, compare->x1); + base->y0 = PS_MIN(base->y0, compare->y0); + base->y1 = PS_MAX(base->y1, compare->y1); + return; +} + +static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System int numCols, int numRows // Size of image ) { psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms - if (!region) { - region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); - } + psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell // Limit the region using the nominated point (X,Y) @@ -25,5 +40,5 @@ fromCoords->x = (X); \ fromCoords->y = (Y); \ - psPlaneTransform(toCoords, wcs->transform, fromCoords); \ + psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \ toCoords->x += wcs->crpix1; \ toCoords->y += wcs->crpix2; \ @@ -33,10 +48,76 @@ region->y1 = PS_MAX(region->y1, toCoords->y); - REGION_RANGE(0,0); // Lower left - REGION_RANGE(0,numRows); // Upper left - REGION_RANGE(numCols,0); // Lower right - REGION_RANGE(numCols,numRows); // Upper right + REGION_RANGE(0, 0); // Lower left + REGION_RANGE(0, numRows); // Upper left + REGION_RANGE(numCols, 0); // Lower right + REGION_RANGE(numCols, numRows); // Upper right return region; +} + +// Activate/deactivate a single element for a list +void fileActivationSingle(pmConfig *config, // Configuration + const char **files, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ) +{ + assert(config); + for (int i = 0; files[i] != NULL; i++) { + pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file + } + return; +} + +// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells +static pmFPAview *filesIterateDown(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + return view; +} + +// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells +static bool filesIterateUp(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + psFree(view); + return true; } @@ -49,8 +130,9 @@ psVector *cdelt1 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT1 values psVector *cdelt2 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT2 values - psArray *regions = psArrayAllocEmpty(BUFFER); // Region for projection + psArray *projRegions = psArrayAllocEmpty(BUFFER); // Region for projection int numProj = 0; // Number of projections psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input + psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image for (int i = 0; i < data->numInputs; i++) { @@ -68,7 +150,4 @@ int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows - double xLow = wcs->crpix1 / fabs(wcs->crpix1) + wcs->cdelt2 / fabs(wcs->crpix2) - double xHigh = numCols + wcs->crpix1 / fabs(wcs->crpix1) + numRows * wcs->cdelt2 / fabs(wcs->crpix2); - pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System if (!wcs) { @@ -76,4 +155,6 @@ return false; } + + psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image bool found = false; // Found a projection? @@ -83,6 +164,7 @@ wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) { - skycellRegion(regions->data[j], wcs, numCols, numRows); + regionMinMax(projRegions->data[j], region); target->data.S32[i] = j; + skycellRegion(wcs, numCols, numRows); found = true; } @@ -94,5 +176,5 @@ psVectorAppend(cdelt1, wcs->cdelt1); psVectorAppend(cdelt2, wcs->cdelt2); - psArrayAdd(regions, regions->n, skycellRegion(NULL, wcs, numCols, numRows)) + psArrayAdd(projRegions, projRegions->n, region); target->data.S32[i] = numProj; @@ -101,2 +183,88 @@ } + pmFPAfileActivate(data->config->files, false, NULL); + + for (int i = 0; i < numProj; i++) { + psRegion *projRegion = projRegions->data[i]; // Region for skycell projection + int xSize = projRegion->x1 - projRegion->x0; // Size of unbinned image + int ySize = projRegion->y1 - projRegion->y0; // Size of unbinned image + // Size of binned image 1 + int numCols1 = xSize / (float)data->bin1 + 0.5, numRows1 = ySize / (float)data->bin1 + 0.5; + // Size of binned image 2 + int numCols2 = numCols1 / (float)data->bin2 + 0.5, numRows2 = numRows2 / (float)data->bin1 + 0.5; + + psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image + psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image + psImageInit(image1, 0); + psImageInit(image2, 0); + + for (int j = 0; j < data->numInputs; j++) { + if (target->data.S32[j] != i) { + continue; + } + + fileActivationSingle(data->config, inFiles, true, j); + pmFPAview *view = filesIterateDown(data->config); // View to readout + if (!view) { + psError(psErrorCodeLast(), false, "Unable to iterate down."); + // XXX Cleanup + return false; + } + + pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j); + pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input + psFree(view); + + pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts + if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) { + psError(psErrorCodeLast(), false, "Unable to rebin image"); + // XXX Cleanup + return false; + } + if (!pmReadoutRebin(bin2RO, bin1RO, data->maskVal, data->bin2, data->bin2)) { + psError(psErrorCodeLast(), false, "Unable to rebin image"); + // XXX Cleanup + return false; + } + + psRegion *imageRegion = imageRegions->data[j]; // Region for image + // Offsets for image on skycell + int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1 + 0.5; + int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1 + 0.5; + int xOffset2 = xOffset1 / (float)data->bin2 + 0.5, yOffset2 = yOffset1 / (float)data->bin2 + 0.5; + + // XXX Completely neglecting rotations + // The skycells are divided up neatly with them all having the same orientation + psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "="); + psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "="); + + psFree(bin1RO); + psFree(bin2RO); + filesIterateUp(data->config); + } + + { + psString filename = NULL; // Filename for image + psStringAppend(&filename, "skycell_%d_1.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, image1, 0, NULL); + psFitsClose(fits); + } + + { + psString filename = NULL; // Filename for image + psStringAppend(&filename, "skycell_%d_2.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, image2, 0, NULL); + psFitsClose(fits); + } + + psFree(image1); + psFree(image2); + } + + // XXX Cleanup + return true; +} Index: trunk/ppSkycell/src/ppSkycellVersion.c =================================================================== --- trunk/ppSkycell/src/ppSkycellVersion.c (revision 23982) +++ trunk/ppSkycell/src/ppSkycellVersion.c (revision 23992) @@ -18,5 +18,4 @@ #include #include -#include #include "ppSkycell.h" @@ -78,5 +77,4 @@ psLibVersionHeader(header); psModulesVersionHeader(header); - ppStatsVersionHeader(header); psString version = ppSkycellVersion(); // Software version