Index: /branches/pap_branch_080207/ppStack/src/ppStackCamera.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackCamera.c (revision 16403) +++ /branches/pap_branch_080207/ppStack/src/ppStackCamera.c (revision 16404) @@ -12,4 +12,5 @@ +#if 0 // Define an output convolved image file static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name @@ -64,5 +65,5 @@ return imageFile; } - +#endif @@ -196,4 +197,5 @@ } +#if 0 // Output convolved files pmFPAfile *outconvImage = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config, @@ -217,5 +219,5 @@ return false; } - +#endif psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, @@ -321,4 +323,5 @@ } + // Output PSF if (havePSFs) { Index: /branches/pap_branch_080207/ppStack/src/ppStackLoop.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16403) +++ /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16404) @@ -16,6 +16,4 @@ // All files in the system static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", - "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT", - "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", @@ -29,13 +27,17 @@ static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 }; -// Input files for the combination -static char *inCombineFiles[] = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 }; - -// Output files for the combination and photometry -static char *outCombineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", - "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", - "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", - "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", - "PSPHOT.INPUT.CMF", 0 }; +// Files required for the convolution +static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 }; + +// Output files for the combination +static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 }; + +// Files for photometry +static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", + "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", + "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", + "PSPHOT.INPUT.CMF", 0 }; + +#define CONVOLVED_ALREADY // Already have the convolution products --- testing @@ -78,4 +80,5 @@ } +#if 0 // Set the data level for a list of files static void fileSetDataLevel(pmConfig *config, // Configuration @@ -98,5 +101,5 @@ return; } - +#endif // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells @@ -150,4 +153,28 @@ return true; } + +#ifndef CONVOLVED_ALREADY +// Write an image to a FITS file +static bool writeImage(const char *name, // Name of image + psMetadata *header, // Header + const psImage *image // Image + ) +{ + assert(name); + assert(image); + + psFits *fits = psFitsOpen(name, "w"); + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open FITS file to write image."); + return false; + } + if (!psFitsWriteImage(fits, header, image, 0, NULL)) { + psError(PS_ERR_IO, false, "Unable to write FITS image."); + return false; + } + psFitsClose(fits); + return true; +} +#endif @@ -263,10 +290,26 @@ } + const char *suffix = "conv.fits"; // Suffix for convolved images; ultimately this will be from recipe + const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root + assert(outName); + psArray *imageNames = psArrayAlloc(num); + psArray *maskNames = psArrayAlloc(num); + psArray *weightNames = psArrayAlloc(num); + for (int i = 0; i < num; i++) { + psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images + psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix); + psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix); + psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix); + imageNames->data[i] = imageName; + maskNames->data[i] = maskName; + weightNames->data[i] = weightName; + } // Generate convolutions and write them to disk psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again + psArray *kernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking for (int i = 0; i < num; i++) { pmFPAfileActivate(config->files, false, NULL); - psArray *files = fileActivationSingle(config, prepareFiles, true, i); + psArray *files = fileActivationSingle(config, convolveFiles, true, i); pmFPAview *view = filesIterateDown(config); if (!view) { @@ -276,7 +319,8 @@ } - pmReadout *readout = pmFPAviewThisReadout(view, files->data[0]); // Input readout + pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout psFree(view); +#ifndef CONVOLVED_ALREADY // Background subtraction, scaling and normalisation is performed automatically by the image matching if (!ppStackMatch(readout, sources, targetPSF, config)) { @@ -287,32 +331,15 @@ } -#ifdef CONVOLUTION_FILES - if (readout->image) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_image.fits", i); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, readout->image, 0, NULL); - psFitsClose(fits); - } - if (readout->mask) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_mask.fits", i); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, readout->mask, 0, NULL); - psFitsClose(fits); - } - if (readout->weight) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_weight.fits", i); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, readout->weight, 0, NULL); - psFitsClose(fits); - } -#endif // CONVOLUTION_FILES + // Write the temporary convolved files + pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image + assert(hdu); + writeImage(imageNames->data[i], hdu->header, readout->image); + writeImage(maskNames->data[i], hdu->header, readout->mask); + writeImage(weightNames->data[i], hdu->header, readout->weight); +#endif cells->data[i] = psMemIncrRefCounter(readout->parent); + kernels->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, readout->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL)); filesIterateUp(config); } @@ -320,15 +347,15 @@ psFree(targetPSF); - // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take - // care of the basic things, and we can just worry about grabbing the readouts by chunks - fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK); - pmFPAfileActivate(config->files, false, NULL); - fileActivation(config, inCombineFiles, true); - fileActivation(config, outCombineFiles, true); - + // Stack the convolved files { + pmFPAfileActivate(config->files, false, NULL); + fileActivation(config, combineFiles, true); + if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) { + fileActivation(config, photFiles, true); + } pmFPAview *view = filesIterateDown(config); if (!view) { psFree(cells); + psFree(kernels); return false; } @@ -341,4 +368,14 @@ } psFree(cells); + + // FITS files + psArray *imageFits = psArrayAlloc(num); + psArray *maskFits = psArrayAlloc(num); + psArray *weightFits = psArrayAlloc(num); + for (int i = 0; i < num; i++) { + imageFits->data[i] = psFitsOpen(imageNames->data[i], "r"); + maskFits->data[i] = psFitsOpen(maskNames->data[i], "r"); + weightFits->data[i] = psFitsOpen(weightNames->data[i], "r"); + } // Read convolutions by chunks @@ -349,19 +386,10 @@ assert(readout); - // Files, containing the FITS handle - pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i); - pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i); - pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i); - - // FITS handles - psFits *fitsImage = imageFile->fits; - psFits *fitsMask = maskFile->fits; - psFits *fitsWeight = weightFile->fits; - - if (!pmReadoutReadChunk(readout, fitsImage, 0, numScans, overlap) || - !pmReadoutReadChunkMask(readout, fitsMask, 0, numScans, overlap) || - !pmReadoutReadChunkWeight(readout, fitsWeight, 0, numScans, overlap)) { + if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap) || + !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) || + !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) { psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i); psFree(readouts); + psFree(kernels); psFree(outRO); psFree(view); @@ -370,40 +398,37 @@ } -#if 0 +#ifdef TESTING + { + pmReadout *ro = readouts->data[0]; + psTrace("ppStack", 1, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols, + ro->row0, ro->row0 + ro->image->numRows); + } +#endif + if (!ppStackReadout(config, outRO, readouts)) { psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); psFree(readouts); + psFree(kernels); psFree(outRO); psFree(view); return false; } -#else - printf("Stack...\n"); -#endif for (int i = 0; i < num && more; i++) { pmReadout *readout = readouts->data[i]; assert(readout); - - // Files, containing the FITS handle - pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i); - pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i); - pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i); - - // FITS handles - psFits *fitsImage = imageFile->fits; - psFits *fitsMask = maskFile->fits; - psFits *fitsWeight = weightFile->fits; - - more &= pmReadoutMore(readout, fitsImage, 0, numScans, overlap); - more &= pmReadoutMoreMask(readout, fitsMask, 0, numScans, overlap); - more &= pmReadoutMoreWeight(readout, fitsWeight, 0, numScans, overlap); - } - } + more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, overlap); + more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, overlap); + more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, overlap); + } + } + psFree(readouts); - - // Get rid of the input files - fileActivation(config, outCombineFiles, false); - filesIterateUp(config); + psFree(kernels); + for (int i = 0; i < num; i++) { + psFitsClose(imageFits->data[i]); + psFitsClose(maskFits->data[i]); + psFitsClose(weightFits->data[i]); + } if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") && @@ -432,5 +457,5 @@ // Write out the output files - fileActivation(config, outCombineFiles, true); + fileActivation(config, combineFiles, true); filesIterateUp(config); Index: /branches/pap_branch_080207/ppStack/src/ppStackReadout.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackReadout.c (revision 16403) +++ /branches/pap_branch_080207/ppStack/src/ppStackReadout.c (revision 16404) @@ -10,4 +10,6 @@ #include "ppStack.h" + +#define ARRAY_BUFFER 16 // Number to add to array at a time #define WCS_TOLERANCE 0.001 // Tolerance for WCS @@ -26,5 +28,5 @@ psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. -// float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution + float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution int num = readouts->n; // Number of inputs @@ -49,13 +51,4 @@ weighting = 1.0; } - -#if 0 - // I don't think 'scale' is needed, since the convolution matches the scales - float scale = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.SCALE"); // Rel. scale - if (!mdok || !isfinite(scale)) { - psWarning("No SCALE supplied for image %d --- set to unity.", i); - scale = 1.0; - } -#endif float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time @@ -135,15 +128,7 @@ #endif - - - // XXX THIS NEEDS WORK -#if 0 for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image pmReadout *readout = data->readout; // Readout for this image - - - // XXX Need some other mechanism to get the subtraction kernel - // Extract the regions and solutions used in the image matching @@ -151,5 +136,5 @@ { psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, - "^SUBTRACTION.REGION$"); // Iterator + PM_SUBTRACTION_ANALYSIS_REGION); // Iterator psMetadataItem *item = NULL;// Item from iteration while ((item = psMetadataGetAndIncrement(iter))) { @@ -159,30 +144,24 @@ psFree(iter); } - psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions + psArray *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels { psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, - "^SUBTRACTION.SOLUTION$"); // Iterator + PM_SUBTRACTION_ANALYSIS_KERNEL); // Iterator psMetadataItem *item = NULL;// Item from iteration while ((item = psMetadataGetAndIncrement(iter))) { assert(item->type == PS_DATA_VECTOR); - solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V); + kernels = psArrayAdd(kernels, ARRAY_BUFFER, item->data.V); } psFree(iter); } - assert(regions->n == solutions->n); - - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, - "SUBTRACTION.KERNEL"); // Kernels - - psPixels *reject = pmStackReject(data->pixels, threshold, regions, - solutions, kernels); // List of pixels to reject + assert(regions->n == kernels->n); + + psPixels *reject = pmStackReject(data->pixels, threshold, regions, kernels); // Pixels to reject psFree(data->pixels); data->pixels = reject; - psFree(solutions); + psFree(kernels); psFree(regions); } -#endif - #ifdef REJECTION_FILES