Changeset 34596 for trunk

Timestamp:

Oct 25, 2012, 11:37:10 AM (14 years ago)

Author:

bills

Message:

New features for ppstamp

-sources <sources file>
-write_cmf

If sources file is supplied (or -write_cmf and astrom is supplied) extract sources
from the region of interest and output them into a cmf file.

Also do a better job insuring that zero point is saved in the output headers.

Location:

trunk/pstamp/src

Files:

• ppstampArguments.c (modified) (4 diffs)
• ppstampMakeStamp.c (modified) (15 diffs)
• ppstampOptions.c (modified) (1 diff)
• ppstampOptions.h (modified) (1 diff)
• ppstampParseCamera.c (modified) (1 diff)
• pstamp.h (modified) (1 diff)

trunk/pstamp/src/ppstampArguments.c

@@ -25 +25 @@
     fprintf(stderr, "   [-mask   mk_image]    : mask image\n");
     fprintf(stderr, "   [-variance var_image] : variance image\n");
+    fprintf(stderr, "   [-sources sources]    : sources cmf (ignored for chip stage)\n");
     fprintf(stderr, "   [-stage stage]        : stage of input image (raw, chip, warp, stack, diff)\n");
     fprintf(stderr, "   [-write_jpeg]         : write a JPEG  format of the image stamp\n");
+    fprintf(stderr, "   [-write_cmf]          : create an output cmf with the sources overlapping the stamp\n");
+    fprintf(stderr, "   [-wholefile]          : ignore the region of interest and process the entire input image\n");
     // fprintf(stderr, "   [-no_censor_masked]   : do not set masked pixels to NAN\n");
     fprintf(stderr, "\n");
@@ -92 +95 @@
         options->stage = psStringCopy(argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
+    } else {
+        // Should we require this? 
+        options->stage = psStringCopy("unknown");
     }
     if ((argnum = psArgumentGet(argc, argv, "-write_jpeg"))) {
         psArgumentRemove(argnum, &argc, argv);
         options->writeJPEG = true;
-    } else {
-        options->writeJPEG = false;
+    }
+    if ((argnum = psArgumentGet(argc, argv, "-write_cmf"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        options->writeCMF = true;
     }
 
-    pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", "-astromlist");
+    // XXX: Note: the various list options have never been tested with ppstamp to my knowledge
+    bool gotAstrom = false;
+    if ((argnum = psArgumentGet(argc, argv, "-astrom"))) {
+        gotAstrom = true;
+        pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", "-astromlist");
+    }
     pmConfigFileSetsMD(config->arguments, &argc, argv, "MASK",   "-mask", "-masklist");
     pmConfigFileSetsMD(config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist");
 
-    // the input file is a required argument; if not found, we will exit
+    if ((argnum = psArgumentGet(argc, argv, "-sources"))) {
+        pmConfigFileSetsMD(config->arguments, &argc, argv, "SOURCES", "-sources", "-sourceslist");
+        // supplying a sources file implies that we want to save the output
+        options->writeCMF = true;
+    } else if (options->writeCMF && !gotAstrom) {
+        // if we didn't get a sources file but the user wanted us to write the sources we must
+        // have an astrometry file supplied
+        // XXX: Is this too restrictive? Could the sources be contained in say the input file?
+        psError(PSTAMP_ERR_ARGUMENTS, true, "cannot write cmf file unless input sources are supplied\n");
+        usage();
+    }
+
+    // the input image file is a required argument; if not found, we will exit
     bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
     if (!status) {
@@ -126 +151 @@
     }
 
-    // finally the output file
+    // finally the only argument left must be output outroot
     if (argc < 2) {
         psError(PSTAMP_ERR_ARGUMENTS, true, "must specify OUTPUT\n");
@@ -135 +160 @@
     }
 
-    // Add the output image (which remains on the command-line) to the arguments list
     psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]);
-
-    // psMetadataPrint(stdout, config->arguments,0);
 
     return config;

trunk/pstamp/src/ppstampMakeStamp.c

@@ -19 +19 @@
 static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
 static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance);
+static bool copySources(pmReadout *outputReadout, pmReadout *inputReadout, pmReadout *astromReadout, psRegion extractRegion);
 
 // convert the input chip's transforms to the output
@@ -80 +81 @@
 }
 
-static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center)
+static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center, pmFPAfile *astrom, pmFPAview *inView)
 {
     pmChip    *outChip;
-    pmFPAview *view = pmFPAviewAlloc(0);
+    pmFPAview *outView = pmFPAviewAlloc(0);
 
     // our output file has a single chip
-    view->chip = 0;
-    outChip = pmFPAviewThisChip(view, output->fpa);
-    psFree(view);
-
-    // copy data from the input chip header to the output.
+    outView->chip = 0;
+    outChip = pmFPAviewThisChip(outView, output->fpa);
+    psFree(outView);
+    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
+
+    // copy data from the input's chip header to the output.
     // since some of the keywords might be duplicated we may not want to copy both
 
-    pmHDU *inHDU  = pmHDUFromChip(inChip);
-    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
+    pmHDU *inHDU  = NULL;
+    if (strcmp(options->stage, "raw") && input != astrom) {
+        // Copy the header from the astrometry file since it contains the psphot and psastro results
+        pmChip *astromChip = pmFPAviewThisChip(inView, astrom->fpa);
+        inHDU  = pmHDUFromChip(astromChip);
+    } else {
+        inHDU  = pmHDUFromChip(inChip);
+    }
 
     if (inHDU->header) {
@@ -102 +110 @@
     }
 
-    // copy the fpa concepts
+    // copy the fpa and chip concepts
     pmConceptsCopyFPA(output->fpa, input->fpa, false, false);
+    // Needed to preserve FPA.ZP
     pmConceptsCopyChip(outChip, inChip, false);
+
+    // copy the cell concepts (needed to get CELL.EXPOSURE which gets used to set EXPTIME)
     pmCell *outCell = outChip->cells->data[0]; // The only output cell
 
@@ -122 +133 @@
     }
 
-    // If input had WCS convert it for stamp
+    // If input had WCS convert it for the new image coordinate system of our stamp
     if (input->fpa->toSky) {
         // copy the input fpa's transforms
@@ -135 +146 @@
         }
     } else {
-        psWarning("No WCS present in input\n");
-    }
-
+        psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "No WCS present in input.");
+    }
+
+    // Save our specific metadata
     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "RA_DEG", PS_META_REPLACE, "Right Ascension of stamp center", RAD_TO_DEG(center->sky->r));
     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "DEC_DEG", PS_META_REPLACE, "Declination of stamp center", RAD_TO_DEG(center->sky->d));
+
+    // XXX: TODO: save the region of interest boundary and the offset to be applied. 
+    // XXX: Gene thinks that there may be a way to handle this in the WCS
+
+    // Now a few sanity checks
+    // Make sure zero point in the fpa concepts matches the observed value in the input header
+    bool status;
+    float zpt_obs = psMetadataLookupF32(&status, outHDU->header, "ZPT_OBS");
+    if (status) {
+        psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.ZP", PS_META_REPLACE, "Magnitude zero point", zpt_obs);
+    }
 
     ppstampVersionMetadata(outHDU->header, options);
@@ -324 +347 @@
 
             if (options->censorMasked && !setMaskedToNAN(config, outReadout->image, outReadout->mask, outReadout->variance)) {
-                 psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
-                 status = false;
-                 break;
+                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
+                status = false;
+                break;
             }
         }
@@ -335 +358 @@
         status = true;
 
+        if (options->writeCMF) {
+            pmReadout *astromReadout = NULL;
+            if (astrom->fpa != input->fpa) {
+                astromReadout = pmFPAviewThisReadout(srcView, astrom->fpa);
+            }
+            if (!copySources(outReadout, astromReadout, astromReadout, extractRegion)) {
+                psError(PS_ERR_UNKNOWN, false, "failed to extract sources from region of interest.\n");
+                status = false;
+                break;
+            }
+        }
+
         psFree(outReadout); // drop reference
     }
@@ -341 +376 @@
 
     if (status) {
-        status = copyMetadata(output, input, inChip, options, center);
+        // For raw stage we need to do more work to build a proper header
+        status = copyMetadata(output, input, inChip, options, center, astrom, view);
     }
     return status ? PS_EXIT_SUCCESS : PS_EXIT_UNKNOWN_ERROR;
@@ -428 +464 @@
     double dy = 0.5 * options->roip.dDEC / fpa->toSky->Ys;
 
-    if (dx > 5000) {
-        fprintf(stderr, "requested width %f too large reducing to 5000\n", dx);
-        dx = 5000;
-    }
-    if (dy > 5000) {
-        fprintf(stderr, "requested height %f too large reducing to 5000\n", dy);
-        dy = 5000;
+    // XXX: why do we limit this?
+    if (dx > 8000) {
+        psWarning( "requested width %f too large reducing to 8000\n", dx);
+        dx = 8000;
+    }
+    if (dy > 8000) {
+        psWarning( "requested height %f too large reducing to 8000\n", dy);
+        dy = 8000;
     }
 
@@ -480 +517 @@
     ppstampOverlap   returnval = PPSTAMP_OFF;
 
-//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));
-
     // set up the astrometry
     pmHDU *hdu = pmFPAviewThisHDU(view, astrom->fpa);
@@ -498 +533 @@
         }
     } else {
-        // we use a default FPA pixel scale of 1.0
         if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0) && mustHaveAstrometry) {
             psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
@@ -537 +571 @@
             center->chip->yErr = 0;
             onChip = true;
-            chipToSky(center, input->fpa, chip);
+            if (input->fpa->toSky) {
+                chipToSky(center, input->fpa, chip);
+            }
         }
     }
@@ -550 +586 @@
                 int width  = options->roip.dX;
                 int height = options->roip.dY;
-                if (width > 7000) {
-                    fprintf(stderr, "requested width %d too large reducing to 7000\n", width);
-                    width = 7000;
+                if (width > 8000) {
+                    psWarning( "requested width %d too large reducing to 8000\n", width);
+                    width = 8000;
                 }
-                if (height > 7000) {
-                    fprintf(stderr, "requested height %d too large reducing to 7000\n", height);
-                    height = 7000;
+                if (height > 8000) {
+                    psWarning( "requested height %d too large reducing to 8000\n", height);
+                    height = 8000;
                 }
 
@@ -675 +711 @@
 
     if (!foundOverlap && (returnval == PS_EXIT_SUCCESS)) {
-        fprintf(stderr, "ROI not found in input\n");
+        psWarning( "ROI not found in input\n");
         returnval = PSTAMP_NO_OVERLAP;
     }
@@ -739 +775 @@
     }
 
-    // XXX: this shouldn't be a fprintf
-    fprintf(stderr, "Censored %ld masked pixels\n", numCensored);
+    psLogMsg("ppstamp", PS_LOG_DETAIL, "Censored %ld masked pixels\n", numCensored);
 
     return true;
 }
+
+static bool copySources(pmReadout *outReadout, pmReadout *inReadout, pmReadout *astromReadout, psRegion extractRegion) {
+    bool status;
+
+    // first look for detections in the input readout. Those would have come from the -sources file
+    pmDetections *inDetections = psMetadataLookupPtr (&status, inReadout->analysis, "PSPHOT.DETECTIONS");
+    if (!inDetections) {
+        // No detections on the input readout. Try the astrometry readout if supplied.
+        if (astromReadout) {
+            inDetections = psMetadataLookupPtr (&status, astromReadout->analysis, "PSPHOT.DETECTIONS");
+            if (!inDetections) {
+                psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no detections found on input or astrometry readout\n");
+                return false;
+            } else {
+                psLogMsg("ppstampMakeStamp", PS_LOG_INFO, "using detections from astrom file\n");
+            }
+        } else {
+            psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no detections found on input readout\n");
+            return false;
+        }
+    }
+    if (!inDetections->allSources) {
+        psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no sources array found on detections structure\n");
+        return false;
+    }
+
+    // Ok we have sources, allocate a detections structure and copy them to the output readout
+    psArray *inSources = inDetections->allSources;
+
+    pmDetections *outDetections = pmDetectionsAlloc();
+    psMetadataAddPtr (outReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN,
+        "psphot detections", outDetections);
+
+    psArray *outSources = outDetections->allSources = psArrayAllocEmpty( 100 );
+
+    // include sources that are off the stamp by just a little bit
+    #define EXPAND_SIZE 16
+    psRegion sourceRegion = extractRegion;
+    sourceRegion.x0 -= EXPAND_SIZE;
+    sourceRegion.y0 -= EXPAND_SIZE;
+    sourceRegion.x1 += EXPAND_SIZE;
+    sourceRegion.y1 += EXPAND_SIZE;
+
+    for (int i = 0; i < inSources->n; i++ ) {
+        pmSource *src = inSources->data[i];
+
+        float x = src->peak->xf;
+        float y = src->peak->yf;
+        if (x < sourceRegion.x0 || x > sourceRegion.x1 ||
+            y < sourceRegion.y0 || y > sourceRegion.y1) {
+            continue;
+        }
+        psArrayAdd(outSources, 100, src);
+    }
+
+    return true;
+}

trunk/pstamp/src/ppstampOptions.c

@@ -32 +32 @@
     options->censorMasked = false;
     options->writeJPEG = false;
+    options->writeCMF = false;
     options->nocompress = false;
     options->wholeFile = false;

trunk/pstamp/src/ppstampOptions.h

@@ -14 +14 @@
     bool        censorMasked;
     bool        writeJPEG;
+    bool        writeCMF;
     bool        nocompress;
     //

trunk/pstamp/src/ppstampParseCamera.c

@@ -120 +120 @@
     }
 
+    // Set up the input and output sources files if needed
+    if (options->writeCMF) {
+        bool status;
+
+        // see if -sources file was supplied. 
+        // If so define the file.
+        psPtr sourcesFile = psMetadataLookupStr(&status, config->arguments, "SOURCES");
+        if (sourcesFile) {
+            pmFPAfile *sources = pmFPAfileBindFromArgs(&status, input, config, "PPSTAMP.INPUT.SOURCES", "SOURCES");
+            if (!status) {
+                psError(psErrorCodeLast(), false, "Failed to load file definition for PPSTAMP.INPUT.SOURCES");
+                return false;
+            }
+            if (!sources && !astrom) {
+                psError(psErrorCodeLast(), false, "Failed to define input sources file");
+                return false;
+            }
+        } else {
+            psLogMsg ("ppstamp", PS_LOG_INFO, "output sources file requested but no input supplied. Will try the supplied astrometry file.\n");
+        }
+        
+        pmFPAfile *outputSources = NULL;
+        outputSources = pmFPAfileDefineSkycell(config, output->fpa, "PPSTAMP.OUTPUT.SOURCES");
+        if (!outputSources) {
+            psError(psErrorCodeLast(), false, "Failed to define output sources file");
+            return false;
+        }
+        outputSources->save = true;
+    }
+
     pmFPAfile *chipImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPSTAMP.CHIP");
     if (!chipImage) {

trunk/pstamp/src/pstamp.h

@@ -26 +26 @@
 
 // values for options mask.
-#define PSTAMP_SELECT_IMAGE      1
-#define PSTAMP_SELECT_MASK       2
-#define PSTAMP_SELECT_WEIGHT     4
-#define PSTAMP_SELECT_CMF        8
-#define PSTAMP_SELECT_PSF       16
-#define PSTAMP_SELECT_BACKMDL   32
-#define PSTAMP_SELECT_INVERSE 1024
-#define PSTAMP_SELECT_UNCONV  2048
-#define PSTAMP_USE_IMFILE_ID 16384
+#define PSTAMP_SELECT_IMAGE         1
+#define PSTAMP_SELECT_MASK          2
+#define PSTAMP_SELECT_WEIGHT        4
+#define PSTAMP_SELECT_SOURCES       8
+#define PSTAMP_SELECT_CMF           8
+#define PSTAMP_SELECT_PSF           16
+#define PSTAMP_SELECT_BACKMDL       32
+#define PSTAMP_SELECT_JPEG          64
+// unused                           128
+// unused                           256
+#define PSTAMP_SELECT_UNCOMPRESSED  512
+#define PSTAMP_SELECT_INVERSE       1024
+#define PSTAMP_SELECT_UNCONV        2048
+#define PSTAMP_RESTORE_BACKGROUND   4096
+// unused                           8192
+#define PSTAMP_USE_IMFILE_ID        16384
 
-#define PSTAMP_NO_WAIT_FOR_UPDATE 32768
+#define PSTAMP_NO_WAIT_FOR_UPDATE   32768
+#define PSTAMP_REQUEST_UNCENSORED  0x10000
+#define PSTAMP_REQUIRE_UNCENSORED  0x20000
 
 #define PSTAMP_CENTER_IN_PIXELS 1