Index: trunk/stac/src/stac.c =================================================================== --- trunk/stac/src/stac.c (revision 3673) +++ trunk/stac/src/stac.c (revision 5743) @@ -18,4 +18,5 @@ // Set trace levels psTraceSetLevel(".",0); + psTraceSetLevel("stac",10); psTraceSetLevel("stac.checkMemory",10); psTraceSetLevel("stac.config",10); @@ -38,5 +39,23 @@ // Read input files - psArray *inputs = stacReadImages(config->inputs); + psArray *headers = NULL; // Array of headers + psArray *inputs = stacReadImages(&headers, config->inputs); + + // Generate masks + psArray *masks = psArrayAlloc(inputs->n); + for (int i = 0; i < inputs->n; i++) { + psImage *image = inputs->data[i]; // Image for which to get mask + psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + if (image->data.F32[y][x] <= config->bad) { + mask->data.U8[y][x] = 1; + } else { + mask->data.U8[y][x] = 0; + } + } + } + masks->data[i] = mask; + } // Read maps @@ -47,5 +66,22 @@ // Get size, if not input if (config->outnx == 0 || config->outny == 0) { - stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps); + stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps); + } + + // Fix up the headers + for (int i = 0; i < headers->n; i++) { + psMetadata *header = headers->data[i]; + psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x", + config->outnx); + psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y", + (int)config->outny); + psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32); +#if 0 + bool mdok = true; // Result of MD lookup + float crpix1 = psMetadataLookupF32(&mdok, header, "CRPIX1"); + if (mdok && !isnan(crpix1)) { + psMetadataAddF32(header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, + "WCS Coordinate reference pixel", crpix1 - +#endif } @@ -58,12 +94,12 @@ // Write error images out to check for (int i = 0; i < inputs->n; i++) { - char errName[MAXCHAR]; // Filename of error image - sprintf(errName,"%s.err",config->inputs->data[i]); - psFits *errorFile = psFitsAlloc(errName); - if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); - } - psTrace("stac", 1, "Error image written to %s\n", errName); - psFree(errorFile); + char errName[MAXCHAR]; // Filename of error image + sprintf(errName,"%s.err",config->inputs->data[i]); + psFits *errorFile = psFitsOpen(errName, "w"); + if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); + } + psTrace("stac", 1, "Error image written to %s\n", errName); + psFitsClose(errorFile); } #endif @@ -72,6 +108,6 @@ psArray *transformed = NULL; psArray *transformedErrors = NULL; - (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL, - config->outnx, config->outny); + (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, NULL, NULL, NULL, + config->outnx, config->outny); psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime); @@ -79,39 +115,55 @@ // Write transformed images out to check for (int i = 0; i < inputs->n; i++) { - char shiftName[MAXCHAR]; // Filename of shift image - char errName[MAXCHAR]; // Filename of error image - sprintf(shiftName,"%s.shift1",config->inputs->data[i]); - sprintf(errName,"%s.shifterr1",config->inputs->data[i]); - psFits *shiftFile = psFitsAlloc(shiftName); - psFits *errFile = psFitsAlloc(errName); - psImage *trans = transformed->data[i]; // Transformed image - psImage *transErr = transformedErrors->data[i]; // Transformed error image - if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); - } - psTrace("stac", 1, "Shifted image written to %s\n", shiftName); - if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); - } - psTrace("stac", 1, "Shifted error image written to %s\n", errName); - psFree(shiftFile); - psFree(errFile); - } -#endif - - - - psVector *scales = NULL; // Relative scales between images - psVector *offsets = NULL; // Offsets between images + char shiftName[MAXCHAR]; // Filename of shift image + char errName[MAXCHAR]; // Filename of error image + sprintf(shiftName,"%s.shift1",config->inputs->data[i]); + sprintf(errName,"%s.shifterr1",config->inputs->data[i]); + psFits *shiftFile = psFitsOpen(shiftName, "w"); + psFits *errFile = psFitsOpen(errName, "w"); + psImage *trans = transformed->data[i]; // Transformed image + psImage *transErr = transformedErrors->data[i]; // Transformed error image + if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); + } + psTrace("stac", 1, "Shifted image written to %s\n", shiftName); + if (!psFitsWriteImage(errFile, NULL, transErr, 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); + } + psTrace("stac", 1, "Shifted error image written to %s\n", errName); + psFitsClose(shiftFile); + psFitsClose(errFile); + } +#endif + + psVector *scales = NULL; // Relative scales between images + psVector *offsets = NULL; // Offsets between images (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff, - config->yMapDiff, config->aper); + config->yMapDiff, config->aper); // Rescale the images (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets); // Set the saturation and bad values psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits - psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits + psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits for (int i = 0; i < transformed->n; i++) { - saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i]; - bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i]; + saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i]; + bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i]; + } + + // Save shifted images + if (config->saveShifts) { + psImage *image = NULL; // Copy of image, to remove NAN for output + for (int i = 0; i < transformed->n; i++) { + char shiftName[MAXCHAR]; // Filename of shifted image + sprintf(shiftName, "%s.shift", config->inputs->data[i]); + psFits *shiftFile = psFitsOpen(shiftName, "w"); + image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32); + (void)psImageClipNaN(image, 0.0); + if (!psFitsWriteImage(shiftFile, headers->data[i], image, 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); + } + psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName); + psFitsClose(shiftFile); + psFree(image); + } } @@ -120,5 +172,5 @@ psImage *combined = NULL; (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated, - bad, config->reject); + bad, config->reject); psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime); @@ -128,24 +180,24 @@ // Write rejection images out to check for (int i = 0; i < rejected->n; i++) { - char rejName[MAXCHAR]; // Filename of rejection image - sprintf(rejName, "%s.shiftrej", config->inputs->data[i]); - - psFits *rejFile = psFitsAlloc(rejName); - if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName); - } - psTrace("stac", 1, "Rejection image written to %s\n", rejName); - psFree(rejFile); + char rejName[MAXCHAR]; // Filename of rejection image + sprintf(rejName, "%s.shiftrej", config->inputs->data[i]); + + psFits *rejFile = psFitsOpen(rejName, "w"); + if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName); + } + psTrace("stac", 1, "Rejection image written to %s\n", rejName); + psFitsClose(rejFile); } // Write out pre-combined image - char preName[MAXCHAR]; // Filename of precombined image + char preName[MAXCHAR]; // Filename of precombined image sprintf(preName, "%s.pre", config->output); - psFits *preFile = psFitsAlloc(preName); - if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", preName); + psFits *preFile = psFitsOpen(preName, "w"); + if (!psFitsWriteImage(preFile, NULL, combined, 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", preName); } psTrace("stac", 1, "Pre-combined image written to %s\n", preName); - psFree(preFile); + psFitsClose(preFile); #endif @@ -153,16 +205,16 @@ psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest for (int i = 0; i < inputs->n; i++) { - regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i], - ((psImage*)(inputs->data[i]))->numCols, - ((psImage*)(inputs->data[i]))->numRows); -#ifdef TESTING - char regionName[MAXCHAR]; // Filename of region image - sprintf(regionName,"%s.region",config->inputs->data[i]); - psFits *regionFile = psFitsAlloc(regionName); - if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName); - } - psTrace("stac", 1, "Region image written to %s\n", regionName); - psFree(regionFile); + regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i], + ((psImage*)(inputs->data[i]))->numCols, + ((psImage*)(inputs->data[i]))->numRows); +#ifdef TESTING + char regionName[MAXCHAR]; // Filename of region image + sprintf(regionName,"%s.region",config->inputs->data[i]); + psFits *regionFile = psFitsOpen(regionName, "w"); + if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName); + } + psTrace("stac", 1, "Region image written to %s\n", regionName); + psFitsClose(regionFile); #endif } @@ -170,23 +222,34 @@ // Transform rejected pixels to source frame psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac, - config->grad, config->inputs); + config->grad, config->inputs); // Get regions of interest in the output frame psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8); for (int y = 0; y < config->outny; y++) { - for (int x = 0; x < config->outnx; x++) { - combineRegion->data.U8[y][x] = 0; - } - } - for (int i = 0; i < inputs->n; i++) { - psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx, - config->outny); - psBinaryOp(combineRegion, combineRegion, "+", region); - psFree(region); + for (int x = 0; x < config->outnx; x++) { + combineRegion->data.U8[y][x] = 0; + } + } + for (int i = 0; i < inputs->n; i++) { + psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx, + config->outny); + psBinaryOp(combineRegion, combineRegion, "+", region); + psFree(region); + + // Do OR of masks + psImage *mask = masks->data[i]; // Mask of input image (bad columns etc) + psImage *cr = rejectedSource->data[i]; // Mask of CRs + for (int y = 0; y < mask->numRows; y++) { + for (int x = 0; x < mask->numCols; x++) { + if (cr->data.U8[y][x] > 0) { + mask->data.U8[y][x] = 1; + } + } + } } // Redo transformation with the masks and scales/offsets - (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource, - combineRegion, scales, offsets, config->outnx, config->outny); + (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, + combineRegion, scales, offsets, config->outnx, config->outny); // Combine the newly-transformed CR-free images, no rejection @@ -194,27 +257,17 @@ rejected = NULL; (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated, - bad, config->reject); + bad, config->reject); // Write output image - psFits *outFile = psFitsAlloc(config->output); - if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); + psFits *outFile = psFitsOpen(config->output, "w"); + int numPix = psImageClipNaN(combined, 0.0); + if (numPix > 0) { + psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix); + } + if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); } psTrace("stac", 1, "Combined image written to %s\n", config->output); - psFree(outFile); - - // Save shifted images - if (config->saveShifts) { - for (int i = 0; i < transformed->n; i++) { - char shiftName[MAXCHAR]; // Filename of shifted image - sprintf(shiftName, "%s.shift", config->inputs->data[i]); - psFits *shiftFile = psFitsAlloc(shiftName); - if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); - } - psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName); - psFree(shiftFile); - } - } + psFitsClose(outFile); // Free everything I've used @@ -228,8 +281,11 @@ psFree(inverseMaps); psFree(errors); + psFree(masks); psFree(transformedErrors); psFree(transformed); psFree(rejectedSource); psFree(combined); + psFree(saturated); + psFree(bad); psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);